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SECTION  I 


INTRODUCTION 

One  aspect  of  aircraft  vulnerability  to  blast  loadings  involves  the 
capability  to  predict  the  dynamic  response  of  aircraft  panels  in  the 
severe  permanent  damage  range.  The  structural  response  program  DEPROP 
was  developed  in  Reference  1 to  calculate  the  linear  elastic  and  elastic- 
plastic,  large  displacement  dynamic  response  of  cylindrical  or  flat 
panels  to  the  pressure  loading  associated  with  a nuclear  blast  wave 
intercepting  an  aircraft  in  flight.  In  this  program  the  loading  was 
assumed  to  be  uniform  over  the  panel  and,  therefore,  only  symmetrical 
response  modes  were  used  in  the  solution.  The  purpose  of  this  study  is 
to  develop  a User's  Manual  for  DEPROP  which  is  applicable  to  general 
non-nuclear  blast  loadings.  This  version  of  DEPROP  computes  strain, 
stress,  displacement  and  reaction  boundary  force  responses  to  an  arbi- 
trary transient  pressure  loading  which  could  excite  both  symmetrical  and 
anti-symmetrical  response  modes.  This  geometric  and  physical  non-linear 
dynamic  response  program  provides  any  combination  of  clamped  and  simply 
supported  boundary  conditions  for  cylindrical  and  flat  panels  which  are 
single  or  multi-layered  with  isotropic  or  orthotropic  material  proper- 
ties. The  DEPROP  program's  elastic  option  can  be  used  with  multiple 
layers  of  isotropic  or  orthotropic  material  while  the  inelastic  option 
is  restricted  to  use  with  only  single  layer  isotropic  panels  based  on  an 
assumed  bilinear  stress-strain  material  representation.  The  DEPROP 
analysis  is  based  on  the  Novozhilov  nonlinear  strain-displacement  rela- 
tions for  large  displacement  response  of  thin  panels  using  the  assump- 
tion of  undeformable  normals.  The  inelastic  formulation  is  based  on  the 
Mises-Hencky  yield  surface,  a kinematic  hardening  model  and  the  Hencky 
stress-strain  relations  from  the  deformation  theory  of  plasticity  with 
modifications  for  regions  of  elastic  unloading  and  reyielding. 

This  User's  Manual  for  the  DEPROP  program  presents  the  theory  and 
program  description  necessary  to  facilitate  the  use  of  the  program. 
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Section  II  presents  the  basic  formulation  of  the  nonlinear  panel  response 
analysis.  The  program  description  and  operation  are  given  in  Section  III 
which  contains  a brief  description  of  each  routine  and  associated  flow 
diagrams,  definition  of  program  variables,  program  input  data  descrip- 
tion and  program  output  description.  The  program  listing  is  presented 
in  Appendix  B. 
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SECTION  II 


PANEL  RESPONSE  ANALYSIS 

The  basic  analytical  formulation  for  DEPROP  is  given  in  Reference  1 
and  is  repeated  herein  for  completeness  of  the  presentation  and  better 
understanding  of  the  DEPROP  program.  The  modifications  and  extensions 
of  DEPROP  involving  the  anti-symmetrical  response  modes,  boundary  reac- 
tion forces,  and  arbitrary  loading  representation  have  been  included  in 
the  presentation  of  the  analytical  formulation. 


2.1  Basic  Theory 

The  single-layered  cylindrical  panel  is  considered  to  have  a 
constant  thickness  h,  mean  radius  a,  subtended  angle  0q  and  length  £. 
The  cylindrical  coordinates  (x,  6,  z)  and  the  axial,  tangential,  and 
radial  displacement  components  (u,  v,  w)  are  shown  in  Figure  1 on  the 
coordinate  surface  which  is  located  at  the  median  surface  of  the  panel. 
The  governing  equations  of  motion  for  the  panel  are  obtained  from  the 
principle  of  virtual  work  for  a dynamic  structural  system  (Reference  2) 
which  is  given  by 


v 


t2  r 


/ /// 

h1  v 


°ij  6eijdV  ' 61  " ff*'  6d  dA 


dt  = 0 


(1) 


where  the  panel  is  undergoing  an  arbitrary  set  of  infinitesimal  virtual 
displacements  6u,  6v,  6w  that  satisfy  the  geometrical  boundary  conditions 
and  vanish  at  t=t^  and  t„;  T is  the  kinetic  energy;  o,  are  the  compo- 
nents  of  total  stress;  are  the  components  of  total  strain;  F is 
the  surface  force  vector;  d is  the  displacement  vector;  and  Integrations 
ar'e  carried  over  volume  V and  deformed  surface  area  A.  It  should  be 
noted  that  this  principle  holds  regardless  of  whether  the  material's 
stress-strain  relations  are  elastic  or  Inelastic  and  whether  the  force 
system  is  conservative  or  nonconservative.  If  it  is  assumed  that 
T = T (u,  v,  w),  then 


6T 


3T 

du 


6u 


. 3T  3T  -. 

+ TT  6v  + rr  6w 
9v  3w 


(3) 
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Figure 


where  the  dots  denote  differentiation  with  respect  to  time.  With  Equation 

(2)  and  using  integration  hy  parts, 


6T  dt  = - 


Vi  ^ r , ViWXx 

77  6u  + — T-r  OV  + 


d 31  . 1 
dt  3*  6W/ 


It  is  assumed  that  the  blast  pressure,  p(x,0,t),  acts  on  the 
coordinate  surface  of  the  cylindrical  or  flat  panel.  As  the  panel 
surface  deforms,  the  elemental  pressure  force  vector  remains  normal  to 
the  coordinate  surface  so  that  it  changes  direction  during  deformation. 
The  magnitude  of  this  force  vector  also  changes  as  the  element  surface 
area  of  the  deformed  panel  changes.  It  should  be  noted  that  the  portion 
of  the  pressure  loading  associated  with  the  force  vector's  dependence  on 
the  deformations  represents  a nonconservative  force  system.  Based  on 
the  rectangular  coordinate  system  (X,Y,Z),  the  components  nY,  nv  and 
nz  of  the  inward  unit  normal  surface  vector  and  the  components  dx,  dy 
and  d^  of  the  displacement  vector  d were  defined  in  Reference  2 in  terms 
of  u,  v and  w and  their  spatial  derivatives.  Thus,  the  vector  dot 
product  of  the  force  and  virtual  displacement  is  expressed  as 


F • 6d  = p(x,e,t)(nxSdx  + nY<Sdy  + nz6dz^ 


By  neglecting  terms  above  the  second  order  and  recasting  in  terms  of  the 
virtual  displacements  6u,  6v  and  6w,  the  virtual  work  done  by  the 
forces  is  given  by 


J J F • 6d  dA  = J /p(x,e,t)(Nu6u  + 6v  + Nw  6w)dA 


where 
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N 

u 

-(w  + w ) 

X X 

N 

V 

-(w0  + wQ  + v)/a 

n 

O 

1 - (w  + w - v.)/a  + u 
8 x 

A 

undeformed  surface  area 

The  subscripts  on  the  displacement  components  denote  spatial  derivatives 

O 

and  w denotes  initial  radial  imperfection  in  the  panel. 

With  Equations  (3)  and  (5)  and  the  relation 


fie 


ij 


Equation  (1)  becomes 
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The  displacement 
series  form  with 


v (t),  w (t); 
mn  mn 


components  are  assumed  in  the  following  truncated 
undetermined  time-dependent  coefficients. 
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u(x,e,t)  = y y u ♦“w*“  (0) 

Z_>  mn  m n 
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v(x,e,t>  = y y v ^(x)^  (e) 
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m=l  n=l 


M N 


w(x,e,t)  = y y » ^(xh* (e) 

f . / mn  m n 


m=l  n*l 


where  <f>  (x)  and  <p  (0)  are  functions  that  satisfy  the  geometric  boundary 
m n 

conditions  of  the  panels.  The  initial  radial  imperfection  in  the  panel 
is  represented  by 

M N 


w(x,0>  = y y a Ox )<t>"( 

/ . / . mn  m n 


m=l  n=l 

where  A are  prescribed  values  based  on  known  or  assumed  deviations 
mn 

from  the  ideal  shape  of  the  panel.  Based  on  Equation  (7),  the  following 
relations  are  obtained: 
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M N 
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6u  - y y Ju  6V  = y y «v  *v , «w  - y y 

/ , / , mn  m n / . / -<  mn  m n f * < -■  mn  m n 


m=l  n*l 


m=l  n*l 


n»“l  n-1 


.u  ,u  3v 
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Introducing  Equation  (9)  into  Equation  (6)  and  since  6umn>  dv^ 

and  6w  are  arbitrary,  the  following  3MN  equations  of  motion  are 
ton 

obtained : 


d_  3T 
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where  the  integrands  of  the  generalized  forces  (2jj^,  2^,  2^)  are  given 
by 


= PN  , 


3u 


pN 


9v 


yw 

’ Smn 


PN 


3w 
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The  kinetic  energy  of  a single-layered  panel  is  given  as 


(ID 


l 0 

T = ^ f J°  (ul 2  + v2  + w2)  dxd6  (12) 

0 0 

where  p is  the  mass  density  of  the  material  and  the  dots  denote  differ- 
entiation with  respect  to  time.  The  rotary  inertia  contributions  to  the 
kinetic  energy  have  been  neglected.  Modification  of  the  mass  density 
for  multilayered  panels  is  introduced  in  subsection  2.5.  Further  development 
of  Equation  (10)  depends  upon  the  establishment  of  the  strain-displacement 
relations,  the  stress-strain  relations,  and  the  displacement  component 
spatial  functions. 
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2.2 


Strain-Displacement  Relations 
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The  strain-displacement  relations  used  in  this  analysis  are 
based  on  the  assumptions:  (1)  strains  are  small  compared  with  unity, 

(2)  the  thickness  of  the  shell  is  small  compared  with  the  radius  and 

(3)  the  Kirchhoff  - Love  hypothesis  that  straight  fibers  which  are 
normal  to  the  undeformed  coordinate  surface  remain  straight  and  normal 
to  the  deformed  coordinate  surface  and  are  not  elongated,  thus  neglec- 
ting transverse  shear  and  normal  strains.  The  basic  formulation  of  the 
following  set  of  nonlinear  strain-displacement  relations  is  attributed 
to  Novozhilov  (Reference  4).  The  total  strain  consists  of  membrane  and 
bending  components  expressed  by  the  form  e = e + zk.  The  membrane 
elongation  and  shear  strains  (e  , eQ_,  e a)  on  the  coordinate  surface 
are  expressed  in  terms  of  the  displacement  components  and  their  spatial 
derivatives: 


1 2 2 2 ° 

e = u +tc[w  +u  +v]+ww 
xx  x 2 x x x xx 


(13a) 


1 1 , . 1 r,  ^ , .2 
Eee  = a Ve  " a Xw  + 7T  t(we  + Xv) 

La 


+ (v0  - Xw)2  + u02]  + w0w0 

a 


(13b) 


ex0  = vx  + I u6  + a wx  (W0  + Xv)  + a vx(v0  " Xw) 


1 1 ° ° 

+ - uQu  + — (w  wQ  + wnw  ) 
a0x  a x0  0 x 


(13c) 


Similarly,  the  change  of  curvature  quantities  (*c  , kqo,  k a)  of  the 

XX  oo  Xu 

coordinate  surface  which  characterize  the  bending  and  torsional  deforma- 
tions of  the  panel  are  given  by 


;| 


i] 

i 

II 

i 

! 


ij 


! 


9 


(14a) 


< * w (1  + vQ/a  - Aw/a  + u ) 

xx  xx  0 x 


*60 = h w90 + h v0  + h (~w  + V + i ux 

a a a 


+ 1 (w00  + Xv  > (v0  " Aw) 


(14b) 


, 1 ..I/ 
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x0  a x0  a x 2 x0  0 x 
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+ WX(W0  + V) 
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(14c) 


Primarily,  only  those  nonlinear  terms  are  included  in  Equation  (14)  which 
involve  the  radial  displacement  and  its  derivatives.  The  subscripts  on  the 
displacement  components  in  Equations  (13)  and  (14)  denote  partial  spatial 
derivatives.  The  end  terms  in  Equation  (13)  are  included  to  account  for 
the  initial  radial  imperfection  of  the  panel  as  indicated  by  Donnell's 
representation  in  Reference  5.  The  parameter  A is  Introduced  in  the 
strain-displacement  relations  so  that  they  apply  to  both  curved  and  flat 
panels.  Thus,  A = 1 for  curved  panels.  For  flat  panels,  A - 0,  a - 1, 

0 is  replaced  by  y,  and  0q  is  replaced  by  b,  the  width  of  the  flat 
panel . 
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2.3  Constitutive  Relations 


In  DEPROP,  the  behavior  of  the  panel  material  is  treated  as 
elastic-plastic  for  isotropic  single-layered  panels  and  elastic  for 
isotropic  and  orthotropic  multilayered  panels.  The  elastic-plastic 
analysis  for  the  single-layered  panel  has  been  established  as  the  basic 
formulation  in  the  DEPROP  program.  The  elastic  multilayered  analysis  is 
established  as  an  alternate  option  based  on  appropriate  modifications  of 
the  elastic-plastic  formulation.  In  the  DEPROP  analysis  the  solution 
involves  total  strains  and  stresses;  therefore,  for  response  in  the 
inelastic  region,  it  is  convenient  to  use  the  deformation  theory  of 
plasticity  instead  of  flow  theory  which  involves  incremental  strains  and 
stresses.  Plastic  deformation  theory  is  based  on  an  averaging  process 
that  permits  a total  strain  solution  dependent  upon  only  the  final 
stress  state  at  the  end  of  a loading  path.  In  general,  deformation 
theory  is  an  approximation  of  the  more  rigorous  flow  (incremental) 
theory  but  is  equivalent  to  flow  theory  for  an  elastic-plastic  material 
when  the  stress  loading  is  proportional,  that  is,  the  ratio  of  principal 
stresses  remain  constant  during  the  loading  process.  However,  since  the 
dynamic  response  solution  is  solved  incrementally  in  time  by  numerical 
methods  in  DEPROP,  the  strain  increments  are  small  and  the  stress  state 
is  fairly  constant  in  the  plastic  region  over  each  time  step  for  which 
the  equations  of  motion  are  solved.  Thus,  the  plastic  deformation 
theory  provides  a much  more  accurate  solution  when  the  averaging  process 
takes  place  separately  over  each  small  time  increment  as  the  response 
solution  is  obtained  by  a step-by-step  timewise  procedure. 

In  deformation  theory  the  total  strain  is  a function  of  the 
state  of  stress  and  consists  of  a recoverable  elastic  component  and  a 
nonrecoverable  plastic  component.  It  is  assumed  that  the  material  is 
incompressible,  that  is,  no  permanent  change  in  volume,  due  to  the 
plastic  strain.  Thus,  the  total  plastic  strain  is  equal  to  the  deviator 
plastic  strain.  Furthermore,  it  is  assumed  that  the  material's  uniaxial 
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stress-strain  curve  is  modeled  by  the  bilinear  representation  shown  in 

Figure  2 in  which  the  strain  hardening  is  defined  by  slope  Et.  This 

stress-strain  representation  is  interpreted  for  the  biaxial  state  of 

stress  through  the  use  of  the  effective  stress  (o)-ef fective  strain  (e) 

concept,  in  which  the  secant  modulus  (E  ) indicated  in  Figure  2 is 

s 

defined  by 


a + E (e  - e ) 
o t o 

e 


(15) 


where  oq,  are  the  yield  stress  and  strain,  respectively,  from  the 

material's  uniaxial  bilinear  representation.  Thus,  the  effective 

stress,  effective  strain*  and  secant  modulus  quantities  are  used  to 

relate  the  biaxial  stress-strain  condition  to  the  assumed  uniaxial 

bilinear  stress-strain  representation  for  the  isotropic  material.  The 

— _ % 

effective  stress  and  strain,  expressed  as  o * f(o^  ) and  e ■ g(e^), 
are  functions  of  the  total  stress  and  strain  components,  respectively, 
and  are  more  conveniently  introduced  in  explicit  form  later  in  the 
development . 

The  Hencky  stress-strain  relations  for  deformation  theory 
(Reference  6)  are  used  in  the  plastic  region  and  are  given  in  the 
following  form: 


Figure  2.  Effective  Stress-Strain  Bilinear  Representation 


into  the  form 
and  are  given  by 


) for  the  case  of  plane  stress  (o  Mo 


where 


v 

s 


ft  (i 

E '2 


v) 


'Xi  1 % a.  'Xj  'Xi  'Xj 

e12  = 2 Gxe’  en  = exx’  e22  = £ee 


It  should  be  noted  that  the  forms  of  the  stress-strain  relations  in  the 
elastic  (Eg=E)  and  plastic  regions  are  identical. 

The  yield  criterion,  in  conjunction  with  a hardening  rule,  and 
the  stress-strain  relation  for  unloading  and  reyielding  by  which  the 
past  strain  history  is  preserved  are  to  be  established.  The  initiation 
of  yielding  for  a biaxial  state  of  stress  is  based  on  the  Mises-Hencky 
yield  criterion  (Reference  6)  given  as 


o 


1/2 


711  °22 


+ 3°12  > 


(18) 


where  a is  the  equivalent  or  effective  stress  and 


11 


xx 


°22  ■ aee*  °12  ~ ax6 


This  yield  criterion  states  that  plastic  flow  will  occur  when  the  equiv- 
alent stress  o reaches  a value  equal  to  the  uniaxial  yield  stress  in 
tension  oq.  A kinematic  hardening  model  is  employed  in  conjunction  with 
the  Mises-Hencky  yield  surface  which  accounts  for  the  Bauschinger  effect 
when  reyielding  occurs  due  to  the  strain  reversals  during  unloading. 

The  Bauschinger  effect  for  a strain  hardening  material  is  described  by 
the  yielding  behavior  of  a material  at  a reduced  yield  stress  when 
reloaded  in  the  opposite  direction  from  that  of  the  initial  yielding. 
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The  kinematic  hardening  models  discussed  in  Reference  7 assume  that 
during  plastic  deformation  the  yield  surface  translates  as  a rigid  body 
in  stress  space  with  the  size,  shape,  and  orientation  of  the  elliptical 
yield  surface  being  invariant.  The  kinematic  hardening  model  to  be  used 
in  this  analysis  is  illustrated  in  Figure  3 for  the  Mises-Hencky  yield 
surface  in  the  plane  of  the  principal  stresses  and  ©2*  Corresponding 
to  the  initial  yielding  position  (i)  and  the  unloading  position  (f) 
indicated  in  Figure  2,  the  rigid  translation  of  the  yield  surface  for  a 
shift  of  the  stress  state  from  position  (i)  to  position  (f)  is  shown  in 
Figure  3.  The  change  in  total  stress  components  from  position  (i)  to 
position  (f)  are  defined  by  and,  similarly,  the  corresponding  change 
in  the  total  strain  components  are  defined  by  ^ , so  that 


% r ^r-l  . r-1  r-1 

“ij  = “ij  + °ij(f)  ‘ °ij(i) 


(19) 


i?  = tsr1 


ij 


ij 


'ur-l 

ij  (f ) " £ij(i> 


where 

r = the  number  of  elastic  unloadings  from  yielded  conditions 
(r=l,2. . .) 

«W 0 

(i)  indicates  initiation  of  yielding  or  reyielding 
(f)  indicates  final  position  prior  to  unloading 

The  yield  criterion  for  the  translated  yield  surface  is  based  on  the 
effective  stress  given  as 


[(On-^i)2  - (.^>(0^)  + (o22^2)2  +3(°12-«J2)2] 


1/2 


(20) 


Furthermore,  it  is  advantageous  in  this  analysis  to  relocate  the  origin 
on  the  e axis  after  each  unloading  such  that  the  extended  elastic 
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* 


unloading  curve  passes  through  the  zero  position.  This  is  accomplished 
by  defining  the  effective  strain  as  follows: 


/■i  2,2 

(1~vs) 


(1-v  +v  ) 

8 S 


- (1-4vs+vs)(^ll 


3 A yr  ^2 
+ (e,o'8io) 


(<'H  - *L>*  + <**-&>*) 

-^l)<'22^2>] 


(21) 


1/2 


d+V 


2 '12  12 


Thus,  the  elastic-plastic  behavior  of  the  material  for  sub- 
sequent yieldings  after  an  unloading  has  occurred  is  always  based  on  the 
same  a versus  e curve  which  originates  at  position  (0,0).  This  approach 
requires  that  the  stress-strain  relations  be  modified  by  the  and 
quantities  for  unloading  and  reyielding  conditions  to  account  for 
the  past  stress-strain  history.  The  general  form  of  the  stress-strain 
relations  for  the  elastic,  elastic-plastic,  elastic  unloading,  and 
plastic  reylelding  regions  are  identical,  so  that  the  general  stress- 
strain  relations  based  on  the  form  of  Equation  (17)  is  given  by 


'vr 

°ij  = “ij  + 


1-v 


[<i-v<V^> + VSk-*Ek>*«] 

(i,j,k  = 1,2) 


(22) 


where,  for  the  following  regions  of  response. 


(a) 

initial  elastic  loading 

(b) 

initial  plastic  loading 

vv  y ■ »„r 

(c) 

qth  elastic  unloading 

v p r 'v  q 

t =*E  , a.  = a ^ 

s ’ ij  “ij 
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(d)  qth  reyielding 


E > 
s 


ij 


q y r 
i.1 ' *ij 


Thus,  there  are  four  basic  regions  of  response  for  which  the  stress- 
strain  relations  have  been  established  by  Equation  (22).  For  an  elastic- 
perfectly  plastic  material,  E^  = 0 and  are  set  equal  to  zero  in 
Equations  (20)  and  (22).  It  should  be  noted  that  for  a strain  hardening 
material,  a stress  path  which  may  move  along  the  yield  surface  (neutral 
loading)  would  not  be  properly  represented  in  the  analysis,  since,  upon 
unloading,  the  yield  surface  would  be  rigidly  translated. 


For  elastic,  isotropic  or  orthotropic  multilayered  panels,  the 
stress-strain  relation  formulation  follows  the  approach  presented  in 
Reference  8.  In  orthotropic  layers,  the  geometric  cylindrical  coordi- 
nate axes  and  principal  orthotropic  direction  are  assumed  parallel.  The 
multilayered  cross  section  for  the  panel  is  shown  in  Figure  4 with  the 
nomenclature  used  in  the  following  formulation.  The  position  of  the 
coordinate  surface  relative  to  the  inner  surface  of  the  panel  is  defined 
by  the  distance  H.  The  membrane  and  bending  stress  resultants  for  the 
multilayered  panel  are  given  by 


°xx  = Cllexx  + C12e00  + FllKxx  + F12K00 
°00  = C22E00  + °12Exx  + F22K00  + F12ICxx 


ax0  ~ C33ex0  + F33Kx0 


*V/  t> 

a - D k + D.-ic.q  + F-.e  + F,„e„. 
xx  11  xx  12  00  11  xx  12  00 


°00  ” D22k80  + D12<xx  + F22e00  + F12£xx 


ax6  = °33Kx0  + F33Ex6 


(23a) 

(23b) 

(23c) 

(23d) 

(23e) 

(23f  ) 


18 


The  elastic  stiffness  constants  and  are  defined 


Bij  (hk  " hk-l) 


■ i 

Fu  ■ I Z)  B1J  [0V  - hk-l>  - 2 H(\  - w 


t IX  [<hk-hk-i>  - 3 r<1y  - \-i> 


+ 35  (\  * Vi> 1 


1-vxv0 


" « ^ 

k k ’ n22 


_!§ - k 

- k k * °33 

1-v  v. 

x e 


G k B k 
x6’  ®12 


. k k , k k 

1-v  v_  1-v  va 

x e x 0 


k k 

E , Eq  are  the  moduli  of  elasticity  in  the  x and  0 directions, 

X u 

respectively,  of  the  kth  layer 

ic  k 

v , v are  Poisson's  ratios  in  the  x and  0 directions,  respec- 

X u 

tively,  of  the  kth  layer 

k 

G is  the  shear  modulus  of  the  kth  layer 

Xu 

h^  is  the  distance  from  the  inner  shell  surface  to  the  outer 

surface  of  the  kth  layer 

n is  the  total  number  of  layers 


For  an  isotropic  material  = EQ  = E,  v^  = vQ  = v and  G^g  «*  2(1 -Pv ) 
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It  has  been  found  that  the  optimal  position  of  the  coordinate 
surface  for  the  most  efficient  modal  convergence  is  at  the  neutral  axis 
of  the  cross  section.  When  the  coordinate  surface  is  located  at  the 
neutral  axis,  the  interaction  stiffnesses  vanish.  These  interaction 
stiffnesses  reflect  the  influence  of  the  change  in  curvature  on  the 
membrane  stress  resultants  and  the  membrane  strains  on  the  bending 
stress  resultants.  As  discussed  in  Reference  8,  for  the  general  case  of 
an  antisymmetrical  orthotropic  multilayered  cross  section,  a neutral 
axis,  which  is  defined  when  all  F^  = 0,  does  not  exist  except  for 
special  combinations  of  elastic  characteristics  of  the  various  layers. 
For  the  general  case,  the  position  of  the  coordinate  axis,  defined  by  H, 
is  established  for  this  analysis  by  setting  F^  = - F^  * F33  * 0 to 

obtain  the  values 


E B«  °v  - hk-?> 

hi  ■ <«> 

2Z  Bu  <ht  - w 

k«i 

and  then  H is  determined  by  averaging  these  values  as  follows: 

H - l (Hu  + + H33)  (25) 

It  should  be  noted  that  for  cases  where  the  neutral  axis  does  exist,  the 
coordinate  surface  is  located  at  this  position  through  the  above  proce- 
dure. When  the  center  of  mass  of  the  cross  section  does  not  coincide 
with  the  neutral  axis,  a slight  discrepancy  in  the  lnplane  inertia  would 
be  introduced  since  the  rotary  inertia  is  not  included  in  this  analysis. 
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2.4  Displacement  Component  Functions 


In  Equation  (7)  the  displacement  components  are  expressed  in 

series  form  as  a product  of  time-dependent  coefficients  and  independent 

spatial  functions  $ (x)  and  <f>  (0).  These  spatial  functions  are  selected 
m n 

so  as  to  satisfy  the  geometric  boundary  conditions  of  the  panels.  The 
boundaries  of  the  panel  are  assumed  to  be  either  clamped  or  simply 
supported  and  spatial  functions  are  defined  to  cover  all  combinations  of 
these  boundary  conditions  for  the  four  edges  of  a panel  defined  by 
x = 0,  i and  0=0,  0q.  On  clamped  edges  the  boundary  conditions 

3w  3w  . 
w = v = u = . = „„  = 0 

3x  30 

are  to  be  satisfied  while  on  simply  supported  edges  the  boundary  condi- 
tions 


are  to  be  satisfied.  Since  the  panels  can  be  non-uniformly  loaded,  the 
assumed  displacement  functions  contain  both  symmetrical  and  anti- 
symmetrical  modes.  The  nondimens ional  variables  y “ and  3 » ^ 
are  introduced  for  use  in  this  analysis.  The  spatial  functions  fo?  the 
u and  v displacements  are  assumed  to  be  the  same  whether  the  edges  are 
clamped  or  simply  supported  and  are  given  by 


(y)  = sin  (m  + l)y 


<f>  (3)  = sin  n3 
n 


<t>  (Y)  = sin  my 
m 


= sin  (n  + 1)3 


r 


. 


The  boundary  combinations  for  the  y and  6 directions  are  based 
on  opposite  edges  being  both  clamped,  both  simply  supported  or  one 
clamped  and  one  simply  supported.  The  w-dlsplacement  functions  for  the 
y and  3 directions  are  based  on  the  natural  vibratory  mode  shapes  of  a 
uniform  beam  with  corresponding  end  boundary  conditions.  These  spatial 
functions  in  the  y and  6 directions  are  given  as  follows  for  the  three 
boundary  combinations: 


For  clamped/clamped  or  clamped/simply  supported  edges 


w 

cosh 

Amy 

m 

Ay 

- cos  

m 

IT 

IT 

w 

cosh 

V 

A 8 

- cos 

n 

TT 

IT 

(27) 


where 


A or  X are  the  roots  of  cos  A.  cosh  A.  = 1 for  the  clamped/ 
m n i i 

clamped  boundary  condition 


A or  A are  the  roots  of  tan  A.  = tanh  A.  for  the  clamped/  simply 
m n i i 

supported  boundary  condition 


cosh  A^  - cos  A^ 


sinh  A^  - sin  A^ 


i = n or  m 


For  opposite  edges  simply  supported 
sin  my 


w 

c 

m 


(28) 


1 


w 

»n 


= sin  n8 


It  should  be  noted  that  the  functions  given  in  Equations  (27)  and  28)  are 
orthogonal. 
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For  the  general  cylindrical  panel,  there  are  nine  combinations 
of  boundary  conditions  provided  by  the  DEPROP  program,  namely. 


1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 

9. 

2.5  Governing  Equations  of  Motion 

With  the  strain-displacement  relations  [ Equations  (13)  and  (14)], 
the  stress-strain  relations  [ Equations  (22)  and  (23)]  and  the  displacement 
component  functions  { Equations  (26)  to  (28)]  defined,  the  governing 
equations  of  motion  ( Equation  (10)]  for  elastic-plastic  deformations 
are  developed  further  by  performing  the  indicated  spatial  integrations. 

For  convenience,  the  dimensionless  quantities  W *=  — , V = — , U=— , 

“ w irx  . ir0  T i.  ir  _ a a a a 

W ■ j , Y “ . p * q~  . L = J = q~  » R “ K = <a  are  introduced 

into  the  formulations.  With  this  notation  and  the  spatial  integration 

for  the  kinetic  energy  in  Equation  (12)  performed  analytically,  the 

governing  equations  of  motion  in  the  radial  displacement  direction 

[ Equation  (10c)]  are  given  by 

I 


Y-direction 
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n 


ir  w 


k k^i2^  + 

Y 6 ran 


^ / / 


h/2 


/ °ij  8W1J  d7dWz  “ / / VYdB 

-h/2 


mn 


(m-l,2,...M) 
(n-1,2, . . .N) 


(29) 


where  for  the  w-equations 

k , k0  /2  for  C-C  or  C-S  opposite  boundaries 

Y P 

= 1//2  for  S-S  opposite  boundaries, 
for  the  u and  v-equations 

k = k ■ 1/ /2~ 

Y 6 

and 

^ - 2L2Rp(l  - XW  - XW  + JV  + J 0 ) fjj— 

1 mn 

\ = -2L2Rp(JVe  + JWp  + XV) 

mn 

$ - -2LRp(W  + W ) 
u Y Y 

mn 


(30) 


Although  Equation  (29)  is  given  In  terms  W , the  equations  of  motion  for 

mn 

the  tangential  and  axial  displacement  directions  can  be  obtained  by 

replacing  W with  V and  U , respectively,  and  using  the  appropriate 
mn  mn  mn 

k , k„,  and  2j  expressions. 

y p 

The  remaining  spatial  integrations  in  Equation  (29)  are  to  be 
accomplished  numerically,  thus  providing  a mechanism  for  discretization 
through  the  spatial  points  selected  to  compute  the  representative  elastic- 
plastic  behavior  throughout  the  panel.  Thus,  a sufficient  number  of 
spatial  points  must  be  specified  to  obtain  a satisfactory  deformation 
response  solution.  For  integration  through  the  thickness  of  the  panel 
in  the  z direction,  it  is  convenient  to  separate  the  integrand 
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into  parts  which  either  are  or  are  not  explicitly  dependent  on  the  z 
variable,  that  is,  involving  membrane  strains  and  bending  strains.  The 
total  strain  quantities  given  in  Equation  (29)  for  an  arbitrary 
position  in  the  panel  consist  of  the  membrane  and  bending  components 
given  by 


ij 


cii + “u 


(31) 


Therefore,  the  integrand  can  be  given  by  fm  + zf^  where 


3e 


xx 


3e 


xx  3W 
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00 


3e 


mm 


00  3W 
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x0 


mn 


x0  3W 


mn 


(32a) 


. 3K 

,b  xx 

f = o 


3K 


xx  3W 


+ o 


00 


3K 


mn 


00  3W 


+ a 


X0 


mn 


x0  3W 


mn 


(32b) 


and  the  total  stress  components  are  obtained  from  Equations  (22)  and 
(19)  in  which  i,j=l  denotes  x and  i,j=2  denotes  0.  The  Legendre-Gauss 
quadrature  formula  (Reference  9)  was  chosen  for  the  numerical  Integration 
in  the  z direction  where  L is  the  number  of  points  selected  through  the 
thickness  of  the  panel.  In  the  y and  0 directions  it  is  convenient  to 
have  even  spacing, and  it  is  advantageous  to  have  spatial  points  on  the 
clamped  edges  and  at  the  center  of  the  panel.  Simpson's  quadrature 
formula  (Reference  9)  possesses  these  desirable  features  and  therefore  was 
selected  over  various  Gaussian  quadrature  formulas.  The  number  of 
spatial  points  selected  in  the  y and  0 directions  are  given  by  M and 
N,  respectively,  where  M and  N must  be  odd  numbers.  By  performing  the 
Indicated  numerical  integrations.  Equation  (29)  is  recast  into  the  form: 
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(m“l,2. . .M) 
(n-l,2...N) 


where 

are  the  zeros  of  the  Legendre  polynomial  P-  (£) 

h r 

zi  " 2 

2(l-cJ) 

Hi  = 2 2 

(L+1)Z 

Hj  or  k “ 4 (j  ’k  ‘ eWn) 

=2  (j,k  * odd,  except  for  j*l,  M or  k * 1,N) 

= 1 (j-1,  M and  k-l,N) 


When  symmetry  is  present  in  both  the  y and  B direction,  for  example, 
only  one-quarter  of  the  panel  need  be  considered  and  the  spatial  inte- 
gration takes  the  form 
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n/2  w/2 

A J J F(y,8)dydB 

o o 


9(M-1)(N-1) 
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I 'jhi 
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2 

S 

k=l 


w(vv 

(34) 


where 

„ . , 4.19  Ell 

Cj  ■ 1 J 1,  2 , « . . 2 


1 iti 

2 J " 2 


For  the  purely  elastic  solution  of  a panel,  the  Integrations  through  the 
thickness  can  be  obtained  analytically,  and  results  In  the  following 
simplification  In  Equation  (33): 


I 

1-1 


Hl£i 


2f 


I 

1=1 


Vi*i 
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(35) 
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where  replaces  o()  in  f [ Equation  (32a)],  o44  replaces  o44  In 


ij 


ij 


ij 


fD  [Equation  (32b)]  and 


m E 

aij  = ~2  L(1“V 

J 1-v 


)eij  + VEkk6lj  . 


b E 


Ij  , 2 

J 1-V 


[a- 


V)Klj  + ^k6lj  J 
(i» j *k  = 1,2) 


(36a) 


(36b) 


The  quantities  e00  and  are  given  by  the  normalized 

versions  of  Equations  (13)  and 
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For  the  elastic  response  solution  of  multilayered  panels,  the 

same  formulation  is  used  with  several  modifications.  The  stiffness 

constants  in  Equation  (23),  C , F and  D , are  considered  to  have  been 
2 3 1J 

divided  by  a,  a and  a , respectively,  and  R * a/h^.  In  Equation  (33), 
p is  replaced  by  p given  by 

f 

n 

? * h~  I pk  (\  - hk-i> 

n k-i 

and 

L 

2 vl  - 2M" 

i-1 

2 wi  ■ 4R2£b 

i-l 

The  a and  a quantities  given  in  Equation  (23)  replace  the  appropriate 
total  stress  quantities  given  in  Equation  (32).  With  these  modifications 
in  DEPROP,  elastic  response  solutions  can  be  obtained  for  multilayered 
panels  of  isotropic  or  orthotropic  material  layers. 
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2.6  Numerical  Analysis 

The  second-order  differential  equations  given  by  Equation  (33) 

and  corresponding  equations  for  V and  U are  to  be  solved  numerically 

mn  mn 

in  time.  The  integration  method  used  to  obtain  an  approximate  timewise 
step-by-step  solution  is  based  on  the  central  difference  formula  given 
by 

Vi  ■ h (4t)2  + 2Xk  - Vi  <3,> 

where 

X represents  the  normalized  undetermined  time-dependent  displace- 
ment coefficients,  W , V and  U 

mn  mn  mn 

At  is  the  time  increment 
k denotes  the  time  step 

In  solving  the  set  of  simultaneous  second-order  differential 
equations,  spatial  integrations  must  be  performed  in  the  y and  g direc- 
tions and  in  the  z direction  for  the  elastic-plastic  solution  during  the 
stepwise  time  integration.  The  required  integrations  are  performed 
numerically  during  each  time  step  using  the  values  of  the  displacement 

coefficients  W , V , and  U for  the  particular  time  step  to  compute 
mn  mn  mn 

the  displacements  and  their  derivatives,  the  strain  quantities,  and  the 
stress  quantities  used  in  Equation  (33) . 
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Several  situations  arise  in  the  implementation  of  the  biaxial 
elastic-plastic  theory  in  DEPROP  which  require  special  numerical  treat- 
ment. These  situations  are  associated  with  the  overshoot  during  the 
time  increment  in  which  yielding  occurs,  the  criteria  for  determining 
elastic  unloading  and  restrictions  if  unloading  is  followed  by  immediate 

reyielding,  and  the  consistent  determination  of  e,  E , v , and  o . during 

s s ij 

each  time  step.  The  special  numerical  schemes  used  to  treat  these  three 
situations  are  described  briefly  in  the  following  three  paragraphs, 
respectively. 


Whenever  a point  in  the  panel  yields  or  reyields  during  a time 
increment  (a  > oq) , the  stresses  cannot,  in  general,  be  computed  on  a 
purely  elastic  basis.  The  computation  of  stresses  should  follow  the 
bilinear  stress-strain  curve,  but  this  is  very  difficult  to  effect  since 


o is  not  a linear  function  of  the  o^'s. 


Instead,  an  iterative  scheme 


is  employed  to  adjust  the  a^'s  proportionately,  so  that  o = oq.  By 
elastic  relations,  the  associated  strains  (o^)  are  determined 

for  the  later  computation  of  $ [Equation  (19)].  The  process  for  cor- 
recting for  overshoot  when  yielding  occurs  between  times  t^  ^ and  t^  is 


illustrated  in  Figure  5.  The  values  of  e and  a are  shown  at  time  t 


k-1’ 


the  values  indicated  at  t^  represent  hypothetical  uncorrected  values. 


By  linearly  adjusting  the  stresses,  the  point  (e 


a ) is  reached, 
o 


It 


is  noted  that  the  actual  point  at  time  t^  should  be  (e^,  o^)  instead  of 
(e^,  0q),  but  the  error  in  stress  is  relatively  small  since  E^/E  « 1. 
The  error  introduced  is  proportional  to  the  size  of  the  integration  time 
step  used. 


For  a point  in  the  panel  which  is  yielding  at  time  t^  elastic 
unloading  is  detected  when,  in  proceeding  to  the  next  time  t^,  the 
equivalent  strain  decreases,  i.e.,  When  this  occurs,  and 

o^  are  recomputed  using  the  elastic  unloading  version  of  Equation (22) . 
Furthermore,  unless  is  less  than  it  is  assumed  that  the  point 

did  not  unload.  This  possible  inconsistency  is  partially  numerical  in 
nature  and  is  partially  due  to  the  nonlinearity  of  the  Equations 


I 


! 


Figure  5.  Correction  for  Overshoot  at  Yielding 


(20)  or  (21).  However,  only  rarely  will  a point  pass  the  strain  criterion 
for  unloading  but  fail  the  stress  criterion.  Due  to  numerical  discrep- 
ancies, it  is  possible  for  the  computed  in  the  unloading  region  to  be 
greater  than  Oq  at  time  t^.  This  inconsistency  results  in  reyielding 
without  any  overshoot  correction  being  made.  In  fact,  this  event  repre- 
sents a numerical  error  and  is  usually  associated  with  the  initial 
stages  of  numerical  instability  of  the  solution.  Consequently,  if  this 
event  frequently  occurs,  the  run  is  automatically  terminated  and  a 
smaller  time  increment  must  be  selected. 


In  the  temporal  integration  sequence,  the  displacement  coefficients 
are  computed  for  the  end  of  the  next  time  step  at  tfcfl  through  the 
central-difference  formula  [Equation  (39)]  given  the  past  displacement 
coefficients  at  t^  and  and  the  acceleration  at  t^.  These  extrap- 

olated displacement  coefficients  are  then  used  to  compute  e . . at  t,  . . 

ij  k+1 

Then,  for  points  in  the  plastic  region,  the  quantities  e,  E , and  v are 

s s 
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evaluated  in  order  to  compute  a at  t^.^.  However,  e is  dependent  on 

v [Equation  (21)],  E is  dependent  on  e [Equation  (15)],  v is  dependent  on 

E [Equation  (17)],  and  a is  dependent  on  E and  v [Equation  (22)]. 

Thus,  a simple  iteration  scheme  is  used  to  simultaneously  solve  for  the 

three  parameters,  e,  E and  v , so  that  a consistent  determination  of 

s s 

can  be  made  at  t^+^.  Then,  with  the  displacements,  their  derivatives, 
and  , all  determined  at  t^+^,  the  accelerations  are  computed  at 

t^^  through  Equation  (33)  and  the  whole  process  is  then  repeated  for  the 

next  time  step. 

A method  has  been  established  for  estimating  time  increments  At  for 
the  temporal  integration  that  would  result  in  stable  solutions  for  the 
majority  of  panel  cases.  The  proper  time  increment  is  a function  of 
geometric  and  physical  properties  of  flat  or  curved  panels  and,  for  the 
DEPROP  formulation,  also  a function  of  the  number  of  modes  used  and  the 
spacing  between  spatial  integration  points.  The  method  for  estimating 
time  increments  is  based  on  formulas  for  the  higher  vibratory  frequen- 
cies of  linear  elastic  panels  which  incorporate  the  aforementioned 
parameters.  The  basic  frequency  formulas  for  single-layered  flat  and 
curved  panels  were  obtained  from  References  10  and  11  and  modified  for 
multilayered  panels.  The  time  increment  is  estimated  by  the  product  of 
the  reciprocal  of  the  frequency  and  an  arbitrary  adjustment  factor. 

Large  displacement  and/or  elastic-plastic  panel  responses  are  nonlinear 
so  that  the  arbitrary  adjustment  factors  are  determined  by  back-figuring 
from  the  time  increments  found  to  give  stable  solutions  for  various 
representative  panels.  In  some  panel  cases  considered  where  the  spacing 
between  integration  points  were  critical  for  numerical  stability,  it  was 
found  that  the  At  formula  used  for  finite  difference  solutions  is  appli- 
cable. This  formula  is  based  on  the  time  for  an  inplane  elastic  wave  to 
propagate  between  mesh  points. 

For  flat  panels  the  governing  time  Increment  used  as  the  initial 
estimate  for  At  in  DEFROP  solutions  is  the  smaller  At  obtained  from  the 
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following  two  formulas.  The  first  formula  is  associated  with  the 
highest  solution  frequency  of  a flat  panel  with  inplane  stresses  which 
are  assumed  to  be  at  a level  corresponding  to  yield  or  ultimate  stress 
of  the  panel  material  and  is  given  by 


^ S.  °max^hk  ^k-l^ 
k-1  


X 

mn 


(40) 


where 


o * yield  or  ultimate  stress 

max 

m * nrt-c 

n * n+c 

and  c*0,  0.15,  and  0.3  for  S-S,  S-C,  and  C-C  boundary  conditions,  respec- 
tively. The  second  formula  is  associated  with  the  elastic  wave  propogation 
between  integration  points  in  the  short  direction  and  is  given  by 


For  curved  cylindrical  panels  the  smallest  At  obtained  from  five 
formulas  is  used  for  the  initial  estimate  for  DEPROP  solutions.  The 
first  two  formulas  are  associated  with  high  frequency  modes  of  cylindri- 
cal shells  and  are  given  by 
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where 


e 
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The  third  formula  is  associated  with  the  elastic  wave  propagation 
between  integration  points  in  the  e-direction  and  is  given  by 
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In  cases  where  the  length  of  the  panel  In  the  x-direction  Is  much 
shorter  than  the  arc  length  in  the  0-direction,  it  was  found  that  the 
flat  panel  formulas  were  more  applicable  than  the  above  three  formulas. 
Thus,  the  fourth  and  fifth  formulas  used  are  based  on  Equations  440)  and 
(41)  with  changing  D22  to  and  b to  a0Q  in  Equation  (40),  and  b to  N 
to  M and  C22  to  in  Equation  (41) . 

The  smallest  time  increment  obtained  from  the  modified  formulas 
used  for  either  flat  or  curved  panels  represents  an  estimated  value  that 
will  generally  give  a stable  solution,  but  does  not  necessarily  repre- 
sent an  optimal  value  for  minimizing  computer  time.  If  a stable  solu- 
tion is  obtained  with  this  estimated  time  increment,  the  time  increment 
can  be  increased  until  two  consecutive  solutions  disagree  significantly 
and  the  penultimate  time  increment  selected  for  future  computer  time 
optimization,  if  desired.  If  the  solution  diverges  using  the  initial 
estimated  time  increment  from  the  formula  procedure,  halving  the  time 
increment  should  easily  result  in  a solution  in  the  stable  range.  The 
DEPROP  program  provides  the  option  of  automatically  using  the  above 
estimated  time  increment  or  having  the  user  select  a value. 

2.7  Boundary  Reaction  Forces 

The  reaction  forces  normal  and  tangent  to  the  coordinate 
surface  along  the  boundaries  of  the  panel  are  determined  in  DEPROP  so 
that  the  loads  on  the  panel  fasteners  can  be  estimated  if  desired.  The 
normal  reaction  forces  per  unit  length  for  a clamped  or  simply  supported 
cylindrical  panel,  Vx  and  VQ,  are  given  in  Reference  12  by 
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where  the  resultant  moments  are  given  by 


h/2  h/2 

M x - / axxZdZ’  M0  “ / a00ZdZ’  Mx0 

-h/2  "h/2 

In  addition,  there  are  concentrated  reaction  forces  acting  at  the 
corners  of  the  simply  supported  panel,  and  they  are  given  by 


h/2 

./  ox0zdz  (46) 
-h/2 


R - 2Mx0. 


(47) 


The  eangential  reaction  forces  per  unit  length  along  the  bound- 
aries are  given  by 
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For  the  elastic  solution  of  the  general  orthotropic  Multi- 
layered panel,  the  reaction  forces  are  computed  in  DEPROP  from  the 
following  expressions: 


V 
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It  should  be  noted  that  for  Equations  (49a)  through  (49c)  the  nonlinear 
terns  of  the  bending  strains  have  been  neglected  in  the  determination  of  the 
normal  reaction  forces. 

In  the  elastic-plastic  solution,  the  inelastic  stresses  are 
complicated  nonlinear  functions  of  strain.  Therefore,  the  partial 
derivatives  indicated  in  Equation  (45)  cannot  readily  be  determined  as 
they  were  for  the  elastic  case.  If  these  derivatives  are  determined 
through  finite  difference  expressions,  however,  the  complexity  of  the 
program  is  not  unduly  increased  and  the  reaction  forces  are  obtained 
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from  stress  quantities  already  computed  in  the  solution  at  the  spatial 
Integration  points  on  and  adjacent  to  the  boundaries  of  the  panel.  Thus, 
the  normal  reaction  forces  for  the  elastic-plastic  solution  are  computed 
in  DEPROP  from  the  following  expressions: 
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Similar  expressions  for  Vx  and  VQ  can  be  defined  for  boundaries  along 
x - o and  e-o,  respectively. 
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The  tangential  reaction  forces  are  given  by 


! X ■, 
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To  estimate  the  loads  on  a fastener  at  a given  location  the 
determined  reaction  forces  per  unit  length  near  the  specified  boundary 
location  can  be  multiplied  by  the  fastener  spacing.  For  the  normal 
reaction  forces  the  sign  convention  designates  positive  force  in  the 
positive  z direction  and  negative  force  in  the  negative  z direction. 

For  the  tangential  reaction  forces  the  sign  convention  designates 
positive  force  in  tension  and  negative  force  in  compression.  It  should 
be  noted  that  the  accuracy  of  the  solution  for  these  reaction  forces 
is  limited  by  the  number  of  modes  and  spatial  integration  points  that 
are  allowed  by  the  core  of  the  computer  used. 


2.8  Approximate  Solution  for  Elastic-Plastic  Response  of 
Sandwich  Panels 

The  elastic-plastic  option  of  DEPROP  is  limited  to  handling 
single-layered  isotropic  panels,  so  that  a thin  three-layered  isotropic 
sandwich  panel  is  reduced  within  the  program  to  an  equivalent  single- 
layered panel  based  on  equating  corresponding  extensional  and  bending 
stiffnesses.  A sandwich  panel  with  face  sheets  of  the  same  material 
described  by  oq,  E and  Et,  can  be  reduced  to  an  equivalent  single- 
layered panel  defined  by  the  following  quantities  in  terms  of  the 
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nomenclature  of  Figure  4 


equivalent  thickness 


elastic  modulus  of  face  sheets  (equal  tq  a /e  ) 


equivalent  elastic  modulus 


strain  hardening  modulus  of  face  sheets 


equivalent  strain  hardening  modulus 


equivalent  mass  density 


equivalent  yield  stress 


2.9 
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Panel  Loading  Program  Options 

Four  transient  pressure  loadings  option  are  incorporated  into 
the  DEPROP  program  to  describe  the  pressure  function  p (x,6,t)  which  is 
designated  as  positive  in  the  negative  z-direction.  The  first  loading 
option  is  an  analytical  representation  of  the  transient  pressure  on  a 
flat  panel  generated  from  the  detonation  of  a projectile  fired  into  an 
adjacent  fluid  medium.  The  second  option  provides  a method  for  a more 
arbitrary  pressure  loading  through  spatial  and  temporal  discretization 
of  the  pressure  field.  The  third  option  contains  a simple  uniform 
pressure  distribution  over  the  panel  with  a combination  of  exponential 
and  triangular  time  decay  behind  a sharp-edged  blast  front.  The  fourth 
option  is  also  for  a uniform  pressure  loading,  but  with  an  arbitrary 
temporal  discretization.  These  loading  options  are  discussed  briefly  in 
the  following  paragraphs. 


2.9.1  Loading  Option  1 

This  option  describes  a transient  pressure  loading  p(x,y,t) 
generated  from  test  data  obtained  on  a flat  plate  subjected  to  the  pres- 
sure loading  from  the  detonation  of  a projectile  fired  into  an  adjacent 
fluid  medium.  As  shown  in  Figure  6 the  trajectory  of  the  projectile  is 
defined  by  the  obliquity  angle  <f>  and  the  detonation  position  is  a distance 
Z from  the  flat  panel.  This  pressure  model  is  based  on  the  assumptions 
that  the  shock  wave  shape  is  spherical,  shock  velocity  is  constant  at 
58,800  inch/second  and  the  pressure  pulse  is  triangular  with  a duration 
of  t^.  If  it  is  assumed  that  time  (sec)  initiates  when  the  blast  wave 
reaches  the  center  of  the  panel,  the  time  of  arrival  of  the  blast  wave 
at  an  arbitrary  point  (x,y)  on  the  panel  is  given  by 


t «^_z 
a 58800 


(53) 
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Figure  6.  Loading  Geometry 


where  R is  the  distance  from  the  detonation  point  given  by 


. 2 . 2 _2. 

(x  + y + Z ) 


1/2 


The  normal  pressure  p(x,y,t)  at  this  arbitrary  point  on  the  panel  is 
described  as: 
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p(x,y,t)  = 0 


<t  < ta> 


p(x,y,t)  =*  0 


<C  " Ca  + td> 


t - t 
a 


p(x,y,t)  = pm  cos  c (1  + — ) (ta  1 t < + td) 

d 


pm  " 484  * 95  (t  + 58800  ) 


(psi) 


cos  a = Z/R 


A - B (t  + 58800  ) 


(sec) 


for  <)»  = 0°  : A - 87  x 10  , B * .0686 

=*  30° : A = 90  x l(f6  , B - .1127 


=60°:  A « 84  x 10  , B - .1275 


2.9.2  Loading  Option  2 

A nonuniform  load  can  be  applied  to  either  a curved  or  flat  panel 
by  specifying  a discrete  pressure-time  history  for  an  array  of  selected 
points  covering  the  surface  of  the  panel.  The  spatial  array  of  points 
need  not  coincide  with  the  integration  grid  point  (determined  by  M,  N) , 
but  must  be  a regular  grid  in  the  sense  that  all  points  remain  in  rows 
and  columns  in  the  x-0  plane,  although  the  spacing  between  rows  (and 
columns)  need  not  be  constant.  The  spatial  grid  should  also  span  the 
entire  portion  of  the  panel  analyzed  or  the  program  will  be  forced  to 
linearly  extrapolate  pressures  toward  the  edges. 


The  timewise  variation  is  specified  at  a set  of  evenly  spaced 
times  - the  spacing.  At,  is  the  same  for  all  points.  However,  the 
time  history  of  each  point  does  not  begin  until  time  corresponding 
to  a unique  delay  time  has  been  reached.  This  delay  time  corresponds 
to  the  time  of  shock  arrival  and  is  specified  on  input  for  each  grid 
point.  It  is  important  that  the  first  point  in  the  array  to  be  engulfed 
have  a delay  time  equal  to  zero. 

Pressures  at  intermediate  times  and  interior  spatial  points 
are  determined  by  linear  interpolation.  No  attempt  has  been  made  to 
estimate  shock  arrival  at  interior  points,  thus,  the  shock  wave  will 
tend  to  be  smeared  unless  a great  number  of  grid  points  are  used. 

For  times  beyond  the  last  time  allowed  for  in  the  loading,  a pressure 
equal  to  the  last  value  specified  at  that  grid  point  is  used. 

2.9.3  Loading  Option  3 

The  third  load  option  assumes  a uniform  distribution  over 
the  surface  of  the  panel,  with  simultaneous  engulfment.  A single 
pressure-time  history  describes  the  entire  loading  sequence.  The 
pressure  loading  is  an  analytical  representation  of  a combination  of 
triangular  and  exponential  decay,  as  indicated  in  Figure  7.  The  pres- 
sures in  the  three  regions  indicated  in  Figure  7 are  given  in  analy- 
tical form  as  follows: 

Pj(t)  “ PX  (1  “ f"  > (t  < t') 

n — 

pII(t)  * Po  (1  ' F~  } e to  (t*  < t < tQ)  (55) 

o 

pm(t>  = 0 (t  i to) 
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Figure  7.  Analytical  Pressure  Time  History 


It  should  be  noted  that  the  second  function  is  used  for  time  greater 
than  or  equal  to  t1;  hence,  by  specifying  t'  ■ 0 the  special  loading 
cases  indicated  in  Table  I (step,  triangular,  impulse,  exponential)  can 
easily  be  generated,  where  I is  the  Impulse  and  At  is  the  integration 
time  interval. 


2.9.4  Loading  Option  4 

Like  the  previous  option,  loading  option  4 is  appropriate  for 
uniformly  applied  loads  without  consideration  of  engulfment.  With  this 
option,  discrete  values  of  pressures  are  specified  at  a set  of  times, 
beginning  at  zero.  For  other  times,  linear  interpolation  is  used, 
except  after  the  last  time  in  the  table,  when  a pressure  equal  to  the 
last  value  is  assumed. 
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SECTION  III 


PROGRAM  DESCRIPTION  AND  OPERATION 


. The  DEPROP  computer  program  determines  the  structural  response  of 
cylindrical  or  flat  panels  to  an  applied  pressure  loading.  The  dynamic 
response  can  be  either  elastic  or  elastic-plastic,  where  the  elastic 
option  applies  to  single  and  multilayered  panels  of  isotropic  or  ortho- 
tropic material,  and  the  elastic-plastic  option  applies  to  single  layer 
panels  of  isotropic  material.  Any  combination  of  clamped  and  simply- 
supported  edge  conditions  is  permitted. 

The  dynamic  pressure  loading  options  provided  in  the  program  are  in 
the  form  of  either  an  analytical  function  of  (x,  y,  t)  or  a discrete 
point-by-point  data  specification,  so  that  either  arbitrary  loadings  or 
certain  specific  loadings  are  accommodated  by  the  program.  Provision  is 
also  made  in  the  program  to  determine  the  elastic  response  to  a uni- 
formly applied  static  load,  or  the  static  load  followed  by  a dynamic 
load. 

The  first  two  subsections  describe  the  code,  including  a brief 
summary  of  individual  routines,  flow  diagrams,  and  definitions  of  the 
major  program  variables.  The  appropriate  dimensions  of  all  program 
variables  are  also  listed  in  case  that  the  user  should  want  to  change 
program  dimensions. 

Subsections  3.3  to  3.5  deal  with  actual  program  operation  - how  to 
input  the  data,  program  requirements,  and  a description  of  the  output, 
including  certain  error  messages  which  might  be  encountered. 

A SUPER  INDEX  listing  of  all  program  variables  is  presented  in 
Appendix  A,  and  a complete  program  listing  can  be  found  in  Appendix  B. 
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3.1  Description  of  the  Routines 


Table  2 lists  the  19  routines  and  all  common  blocks  which  make  up 
DGPROP.  The  decimal  length  of  each  common  is  also  Indicated. 

Flow  diagrams  of  the  major  routines  are  presented  In  Figures  8 
to  12,  while  brief  descriptions  of  the  purpose  of  all  the  routines  are 
given  below.  Also  included  are  lists  of  the  routines  which  are  refer- 
enced by  other  routines  and  vice-versa. 


DEPROP 

Main  program.  Reads  preliminary  input  data  and  controls  program  flow. 
Calls  PINIT,  PROP. 

BOLT 

Sets  up  W mode  shapes  for  boundary  conditions  selected. 

Called  by  DSET1 . 

DERV2 


Computes  strains,  displacements,  and  accelerations  in  the  main  inte- 
gration loop. 

Calls  LIST1,  LIST2,  PRESS,  REIT,  SIGMA. 

Called  by  PROP. 

DSET1 


Reads  DEPROP  Input  data  and  calculates  constants. 
Called  by  PROP. 

DSET2 


Calculates  constants  used  in  DEPROP. 
Calls  LEGEND,  DTSTEP. 

Called  by  PROP. 

DSET3 


Calculates  additional  constants  and  writes  out  a description  of  input 
data. 

Calls  BOLT. 

Called  by  PROP. 

DTSTEP 


Computes  an  integration  time  step  small  enough  to  avoid  numerical 
instabilities. 

Called  by  DSET2 . 
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TABLE  2.  DEPROP  ROUTINES  AND  COMMON  BLOCKS 
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E 
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X 
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X 
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X 
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X 

X 

X 
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1 
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X 

X 

X 

X 

X 
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0V013 

1072 

X 

X 

X 

VA0N3 

142 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

ISXIi 

X 

X 

COMMON 

Length 

of 

Common 

Routine 

DEPROP 

BOLT 

DERV2 

DSET1 

DSET2 

DSET3 

DTSTEP 

HIM 

LEGEND 

LIST1 

LIST2 

PINIT 

1 PRESS 

1 PROP 

REIT 

RELAXP 

SEC 

VHDIS  j 

SOLVE 

I 
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START 


IMPUT  NCASES 


Figure  9.  Program  DEPROP  Flow  Diagram 


i 
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(INITIALIZE ) 


(STATIC) 


(DYNAMIC) 


Figure  10.  Subroutine  PROP  Flow  Diagram 
56 
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Figure  11.  Subroutine  DERV2  Flov  Diagram 
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Figure  12.  Subroutine  SIGMA  Flow  Diagram 
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HIM 


! 


Numerical  timewise  integration  routine. 
Called  by  PROP. 

LEGEND 


Sets  up  constants  for  Gaussian  Integration  through  the  thickness  for  an 
elastic-plastic  solution. 

Called  by  DSET2. 

LIST1 

Output  routine  for  the  elastic-only  option. 

Called  by  PROP,  DERV2. 

LIST2 

Output  routine  for  the  elastic-plastic  option. 

Called  by  PROP,  DERV2. 

PINIT 


Reads  in  pressure  loading  time  history.  Sets  up  loading  fjnctions. 
Called  by  PROP. 

PRESS 

Calculates  pressure  loading  at  a given  time. 

Called  by  DERV2. 

PROP 

Executes  the  static  and  dynamic  panel  solutions. 

Calls  DERV2,  DSET1,  DSET2,  DSET3,  HIM,  LIST1,  LIST2,  RELAXP,  SEC. 
Called  by  DEPROP. 

REIT 


Computes  reaction  forces  at  the  boundaries  and  comers. 
Called  by  DERV2. 

RELAXP 


Solves  simultaneous  nonlinear  equations  representing  preblast  condi- 
tions using  a relaxation  procedure. 

Calls  SOLVE. 

Called  by  DEPROP. 
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SEC 

Finds  elapsed  CP  time. 

Calls  system  routine  SECOND. 

Called  by  PROP. 

SIGMA 

Computes  stresses  in  the  main  integration  loop  for  the  elastic-plastic 
option.  Is  not  needed  for  the  elastic-only  option. 

Called  by  DERV2 . 

SOLVE 

Solves  a set  of  simultaneous  linear  algebraic  equations. 

Called  by  RELAXP. 

3.2  Definition  of  Major  Program  Variables 

The  major  program  variables  are  defined  in  this  subsection.  These 
variables  are  listed  alphabetically  with  a brief  description  devoted  to 

each  one.  An  asterisk  preced.  ng  a variable  name  indicates  that  the 

variable  is  input  as  run  data.  The  dimension  of  a variable  is  given 
parenthetically  after  the  variable  name,  where  a numerical  dimension 
indicates  the  fixed  amount  of  storage  required  for  the  variable.  There 
is  no  need  to  change  the  dimension  of  such  a variable.  However,  other 
variables  have  variable  dimensions.  As  an  example,  one  of  the  first 
dimensioned  variables  listed  is  BETR(NBAR) . The  dimension  is  a variable, 
NBAR,  which  represents  the  number  of  spatial  integration  points  used  in 
the  beta-direction  in  the  panel.  This  dimension  must  be  the  largest 
number  of  points  the  user  intends  to  employ.  In  the  DEPROP  program, 
this  dimension  is  23.  If  additional  integration  points  are  required, 
the  dimension  of  NBAR  must  be  increased,  thus  increasing  the  dimensions 
for  all  variables  with  the  dimension,  NBAR.  The  current  dimensions 
provided  for  the  variables  in  the  DEPROP  routines  are  given  in  Table  3. 

Almost  all  of  the  variables  which  may  require  dimension  changes  as 
indicated  above  are  contained  in  the  COMMON  blocks.  There  are  a few 
exceptions  and,  in  such  cases,  the  subroutine  in  which  the  variable  is 
dimensioned  is  indicated  in  the  list  of  variables  or  in  Table  4.  If 
the  dimensions  are  changed,  certain  additional  change-  in  the  program 
may  be  required.  These  changes  are  also  indicated  in  Table  4. 
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Many  of  the  variables  found  in  the  program  result  from  their  use  in 
the  larger  NOVA-2  code  (Reference  1).  In  most  cases,  these  variables 
have  little  or  no  use  in  DEPROP  and  are  so  indicated. 
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TABLE  3.  DIMENSIONS  OF  DEPROP  VARIABLES 


VARIABLE 

DIMENSION 

LBAR 

6 

MB 

13 

MBAR 

23 

MBAR*  NBAR*LBAR1 

1805 

MG 

13 

MGMB2 

49 

NBAR 

23 

NGNBT3 

361 

NKP 

46 

NL 

8 

NPX 

10 

NPY 

10 

1.  This  constraint  is  only  significant  for  an  elastic-plastic  run 

(NDERV=2).  Three  possible  combinations  using  maximum  dimensions 
are:  (17x17x6),  (19x19x5),  (21x21x4). 

2.  The  total  number  of  modes  selected  (MGMB)  from  the  total  possible 
(MB*MG)  cannot  exceed  49. 

3.  The  total  number  of  spatial  Integration  point  allowed  (NGNBT) 
from  the  total  possible  (MBAR*NBAR)  is  361. 


TABLE  4.  PROGRAM  CHANGES  REQUIRED  BY  DIMENSION  CHANGES 


When  Changing  the 
Dimensions  Corre- 

Also  Change  the  Fixed-Point  Number  in  the 
Indicated  Statement 

sponding  to: 

Subroutine 

Location^" 

MGMB*3 

HIM 

COMMON  BLOCK 

LBAR 

LEGEND 

300-11 

MGMB*3 

RELAXP 

20+1 

The  location  code  is  read  as  follows:  s refers  to  the  nth  line 

after  tatement  number  s. 


Major  Program  Variables  for  DEPROP  Routine 


*A 


*AA 

ALTT(LMAX) 

ALXT(LMAX) 

ALXX(LMAX) 

*ANN 


AZ 

BE1 (LMAX) 
BE2 (LMAX) 
BE3(LMAX) 
BETR(NBAR) 
BTL(NL) 

BXL(NL) 

CCRIT(NL) 

CCl(MG) 

CC2 (MB) 

CC5 (MG) 
CC6(MB) 
CINST(3) 
CK(6) 

CM11,CM12 , 
CM2  2,  CM3  3 


Radius  of  the  cylindrical  panel,  a,  in.  Is 
set  equal  to  1.0  for  flat  panel. 

Pressure  loading  parameter,  a,  dimensionless 
(Load  Option  3) 

'V/ 

Change  in  stress  component,  a , psi 

C7U 

Change  in  stress  component,  a^g,  psi 

Change  in  stress  component,  a , psi 

Pressure  loading  parameter,  n,  dimensionless 
(Load  Option  3) 

Constant  equal  to  a/tQ,  1/sec  (Load  Option  3) 
Change  in  strain  component,  $ , in/in. 

Change  in  strain  component,  , in/ in. 

OO 

Change  in  strain  component,  $xg»  in/ in. 

Integration  positions  in  the  beta-direction,  rad 

k.  k k 

Constants  used  in  elastic  option,  E„/(l-v  v„), 

. 9x0 

psi 

k k k 

Constants  used  in  elastic  option,  E / (1-v  v_), 

. x x 0 

psi 

Not  used 

Constant , m,  dimensionless 
Constant,  n,  dimensionless 
Constant,  m+1,  dimensionless 
Constant,  n+1,  dimensionless 
Not  used 

Constants,  1/k  , 1/k  , for  the  equations  of 
motion,  dimensionless 

Stiffness  constants  C^  used  in  elastic,  multi- 
layer integration  through  the  thickness,  lb/in. 
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CN10 

CN11 

CN6 


2 2 
Constant,  2L  R for  elastic  option,  2L  for 

elastic-plastic  option,  dimensionless 

2 2 2 
Constant,  2L  R for  elastic  option,  L /6R  for 

elastic-plastic  option,  dimensionless 

2 

Constant  for  elastic-plastic  option,  E/(l-v  ), 
psi 


CN7 

CN8 

CN9 

COSB (MB*NBAR) 

COSG (MG*MBAR) 

C0S2B (MB*NBAR) 

C0S2G (MG*MBAR) 

CRIT(3) 

DC 

DELT 

*DELTIM 

DELX 

DELX(3*MGMB) 

*DET(NPY,NPY) 

DM11,  DM12, 
DM22,  DM33 

DPRT 

DPRT1 


Shear  modulus  for  elastic-plastic  option, 
E/2 (1+v) , psi 

2 

Constant , L , dimensionless 
2 

Constant,  L /2R,  dimensionless 

Cosine  functions  of  the  beta  angles, 
cos  ((n+l)Bj),  dimensionless. 

Cosine  functions  of  the  gamma  angles, 
cos  ((nri-l)Y^) , dimensionless 

Cosine  functions  of  the  beta  angles, 
cos  (n  B^  ) , dimensionless 

Cosine  functions  of  the  gamma  angles, 
cos  (m  y^) , dimensionless 

Not  used 


Not  used 

Grid  point  spacing  in  B-direction,  Ay,  in., 
or  aA0 , in-rad 

Integration  time  interval,  sec 

Grid  point  spacing  in  Y-direction,  Ax,  in. 

Working  array  in  RELAXP,  dimensionless 

Delay  time  to  when  spatial  point  is  first 
engulfed  by  pressure  wave,  sec  (Load 
Option  2) 

Stiffness  constants  D^  used  in  elastic, 

multilayer  integration  through  the  thickness, 
lb-in. 

Running  time,  in  multiples  of  DPRT1,  used  to 
flag  next  printout,  sec 

Time  interval  between  printouts,  sec 
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*DTIM 

DWB(NGNBT) 

DWG (NGNBT ) 

DWO (NGNBT) 

DX1 (MBAR) 

DY1 (NEAR) 

EC 

EL 

*EM(NL) 

ENX(2*NBAR-2) 

ENT(2*MBAR-2) 

*EP 

EPB(LMAX) 

EPBO(LMAX) 

EPO 

*EPSIF 

ERR(3*MGMB) 


Time  interval  between  specified  pressure 
data,  sec  (Load  Option  2) 


Values  for  imperfection-related  partial 

O 

derivative  W , dimensionless 

p 


Values  for  imperfection-related  partial 

O 

derivatives  W , dimensionless 


Values  for  imperfection-related  displacement 

O 

W,  dimensionless 


Fractional  distance  in  y-direction,  locating 
grid  point  between  two  pressure-mesh  points, 
dimensionless  (Load  Option  2) 

Same  as  DX1,  only  in  B-direction 


Not  used 


Modulus  of  elasticity  for  elastic-plastic 
option,  E,  psi 

Modulus  of  elasticity  for  each  layer  for 
elastic-plastic  option,  E , psi 

Tangential  reaction  force  per  unit  length 
along  boundary,  N^,  lb/in. 

Tangential  reaction  force  per  unit  length 
along  boundary,  N^,  lb/ in. 

Strain  hardening  slope  for  elastic-plastic 
option,  E^,  psi 

Temporary  storage  used  for  either  e or  o 

Equivalent  stress,  squared,  when  response  is 

still  elastic  for  the  elastic-plastic  option, 

-2  ,,2,,  4 

o , lb  /in  . 

Equivalent  yield  strain  corresponding  to 
SIGO  for  elastic-plastic  option,  cq,  in/ in. 

Ultimate  (fracture)  strain  for  elastic-plastic 
option,  e^,  in/in;  not  used 

Allowable  error  in  displacement  coefficients 
in  the  static  solution 


i 


A 
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*ET(NL) 


ETTl(LMAX) 

*EX(NL) 


EXT 

EXT1 (LMAX) 
EXX 

EXXl(LMAX) 

*FG(MB,MG) 

FM11 , FM12 
FM22,  FM33 

FP1 (MG*MBAR) 
FP2 (MG*MBAR) 

FP3 (MG*MBAR) 

FP4 (MG, 2) 

FP5 (MB*NBAR) 
FP6 (MB*NBAR) 
FP7 (MB*NBAR) 


Modulus  of  elasticity  in  the  theta-direction 

v 

for  elastic  option,  E„,  psi 

Temporary  value  of  strain,  e“  , in/ in. 

00 

Strain  component  at  the  time  of  last  yielding, 
e00 * in/in- 

Modulus  of  elasticity  in  x-direction  for  elas- 
tic  option,  E^,  psi 

Temporary  value  of  strain,  em  , in/ in. 

X0 

Strain  component  at  the  time  of  last  yielding, 
ex0,  in/in. 

Temporary  value  of  strain,  e^,  in/ in. 

Strain  component  at  the  time  of  last  yielding, 

e , in/in. 
xx 

Modal  displacement  coefficients  for  the  initial 
imperfections,  in. 


Stiffness  constants  F^  used  in  elastic, 

multilayer  integration  through  the  thickness, 
lb 


Displacement  function,  (x),  dimensionless 

m 

d W 

t—  f (x) , dimensionless 
oj  in 


d2  W 

<t>  (x)  , dimensionless 

3y2  m 


3^m 

at  Y = 0,  tt,  dimensionless 

3Y3 

Displacement  function,  $w  (0),  dimensionless 

n 

3 <j>"(0) 

— gjj , dimensionless 

32^(0) 

, dimensionless 

3B2 
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FP8(MB,2) 

GAM(MBAR) 

GAMMA(41) 

GC 

GX(LBAR) 

*GXT (NL) 

H 

HBAR 

HGO(LBAR) 

*HM(NL) 

ICOMP 

ICOUNT 

IFIRST 

INOUT 

I1IZ(2) 

IP (3*MGMB) 
IXI(MSAR) 

JFIRST 

JL 


a3^(e) 

at  8 = 0,  v,  dimensionless 

3B3 

Integration  positions  in  the  y-direction,  rad 
Not  used 
Not  used 

Zeroes  of  the  Legendre  polynomial  for  the 
elastic-plastic  Gaussian  integration  through 
the  thickness,  dimensionless 

Shear  modulus,  G psi 
x0  r 

Thickness  of  cross  section  for  elastic-plastic 
option,  in. 

Distance  from  the  inner  panel  surface  to  the 
coordinate  surface,  H,  in. 

Weighting  factors  for  the  elastic-plastic 
Gaussian  integration  through  the  thickness, 

H^,  dimensionless 

Distance  from  the  inner  panel  surface  to  the 
outer  surface  of  layer,  h,  in. 

Not  used 

Counter,  initially  set  at  37777000000000000000B, 
and  incremented  for  each  integration  step  in 
program 

Code  designating  first  pass  through  subroutine 
SIGMA,  dimensionless 

Not  used 

Code  designating  the  appropriate  layer  number 
corresponding  to  the  two  panel  surfaces,  dimen- 
sionless 

Working  array  in  RELAXP,  dimensionless 

Integer  locating  grid  points  in  y-direction 
relative  to  pressure-mesh  (Load  Option  2) 

Code  which  indicates  if  the  panel  has  yielded 
for  the  elastic-plastic  option 

Lower  index  on  timewise  interpolation  table 
(Load  Option  4) 
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JLT(NPY.NPX) 
JYJ (NBAR) 

•> 

KALT 

KB 

KC 

KDAM 

*KDS 

KERR 

KOK 

*KPG(NKP) 

*KPB (NKP) 
KSUMA(NGNBT) 

*KTYPE 

KY(LMAX) 

KZ 

*LBAR 

UMAX 

*MB 


Lower  Index  on  timewise  Interpolation  table 
(Load  Option  2) 

Integer  locating  grid  points  in  8-direction 
relative  to  pressure-mesh  (Load  Option  2) 

Not  used 

Not  used 

Not  used 

Not  used 

Response  option  code 

Dynamic  response  error  code  - 0,  no  error; 

1,  error 

Not  used 

Mesh  point  number  (y),  when  paired  with  KPB, 
calling  for  printout 

Mesh  point  number  (8),  when  paired  with  KPG, 
calling  for  printout 

Number  of  z points  which  have  not  yielded  at 
each  spatial  station;  used  only  in  elastic- 
plastic  option 

Code  designating  panel  type 

1,  single-layer  metal 

2,  honeycomb  panel  (3  layers) 

5,  multilayer  panel  (elastic  response  only) 

Code  in  elastic-plastic  response,  indicating 
number  of  times  an  integration  point  has 
yielded,  unloaded,  etc. 

Code  deciding  whether  the  output  routine 
should  print 

Number  of  integration  points  through  the 
thickness ; is  assumed  to  be  one  for  elastic 
option 


Maximum  number  of  integration  points  being 
used ; equal  to  MBAR*NBAR*LBAR 

Number  of  beta  modes  to  be  incorporated  into 
the  solution 
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*MBAR 

*MG 

*MGM(MG) 

MGMB 

MGMB2 

*MOUT(MG*MB) 
MUSE (MB, MG) 
NB 

*NBAR 

*NBN (MB) 
*NBND 
NBT 
NCALL 

NCASE 

*NCASES 

NCHPT 

*NDBUG 

*NDERV 


Number  of  beta  integration  points  to  be  used; 
for  a symmetric  boundary  condition,  only 
approximate  half  as  many  points  are  required 

Number  of  gamma  modes  to  be  incorporated  into 
the  solution 

Gamma  mode  numbers 

Constant,  equal  to  the  total  number  of  modal 
combinations  used,  MG*MB-NNOUT 

Constant,  2*MGMB 

Gamma  modes  not  to  be  included,  in  combination 
with  the  corresponding  NOUT  mode 

Code  designating  which  modal  combinations  are 
to  be  used 

Beta-point  number  corresponding  to  centerline 
of  panel 

Number  of  gamma  integration  points  to  be  used; 
for  a symmetric  boundary  condition,  only 
approximately  half  as  many  points  are  required 
as  otherwise. 

Beta  mode  numbers 

Boundary  condition  code 

Total  number  of  beta  points  monitored 

Code  describing  program  phase: 

0,  find  dynamic  response 

1,  find  static  solution 

2,  read  data,  set  up  constants 

Case  number  currently  being  run 
Number  of  cases  to  be  run 
Not  used 

Output  debugging  control: 

0,  no  debug  output  (normal  option) 

1,  most  debug  output 

2,  all  debug  output. 

Response  code:  1-elastic  (multilayer,  ortho- 

tropic, 2-elastic-plastic  (single  layer, 
isotropic) 
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NELP 

Response  code  for  the  elastic-plastic  option: 
1-keep  solution  elastic,  2-allow  solution  to 
go  elastic-plastic 

NG 

Gamma-point  number  corresponding  to  center 
line  of  panel 

NGNB 

The  spatial  integration  point  number  corre- 
sponding to  the  center  of  the  panel 

NGNBT 

Total  number  of  spatial  integration  points 
used 

NGT 

Total  number  of  gamma  points  monitored 

*NKP 

*NL 

*NLOAD 


NLZ(2*NL) 

NMASS 

*NNOUT 

*NOUT(MG*MB) 

*NPLT 

*NPX 

*NPY 

NREG 

*NSYMB 


Number  of  spatial  points  for  which  printout 
of  strains,  stresses,  reactive  forces,  dis- 
placements and  pressure  is  required 

Number  of  layers . 

Pressure  load  option  code: 

1,  nonuniform,  functions 

2,  nonuniform,  discrete 

3,  uniform,  functions 

4,  uniform,  discrete 

Layer  number  corresponding  to  each  layer's 
upper  and  lower  surfaces 

Not  used 

Number  of  modal  combinations  (<MG*MB)  to  be 
eliminated  from  the  solution 

Beta  modes  not  to  be  included,  in  combination 
with  the  corresponding  MOUT  mode 

Panel  type  code:  0-flat,  1-curved 

Number  of  mesh  points  in  y-directions  for 
which  pressure  data  is  provided  (Load 
Option  2) 

Number  of  mesh  points  in  3-direction  for 
which  pressure  data  is  provided  (Load 
Option  2) 

Not  used 

Symmetry  code  in  the  beta-direction: 
O-symmetric,  1-not  symmetric 
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*NSYMG 


Same  as  NSYMB,  except  in  the  gamma-direction 


NTECO 
* NT  IKE 

NTRIAL 

NU 

NUSE(NBAR.MBAR) 

NY  2 
NZP 
OTTO 

OTT1 

P(NGNBT) 

PB(40) 

PDAM 

*PHI 

PI 

PIMA(MBAR) 

PINA(NBAR) 

PPP 

*PP1 

*PPO 

*PRINT 


Not  used 

Number  of  times  for  which  pressure  data  is 
provided.  (Load  Options  2 and  4) 

Not  used 

Code  indicating  whether  loading  is  spatially 
uniform  or  not:  0-not  uniform,  1-uniform 

Use  code  for  the  spatial  integration  stations 
0-not  used,  1-not  used  for  printout  only, 
2-used  for  integration  only,  3-both 

Constant , equal  to  3*MGMB 

Not  used 

Constant  equal  to  1/t  , 1/sec  (Load 
Option  3)  ° 

Constant  equal  to  1/t  , 1/sec  (Load 
Option  3) 

Pressure  at  each  spatial  point,  psi 
Not  used 
Not  used 

Angle  projectile  trajectory  makes  with  the 
normal  to  the  panel  (z-axis),  <j>,  degrees 
(Load  Option  1) 

Constant,  equal  to  it 

Constants  associated  with  the  equation  of 
motion  and  Simpson's  Rule;  Gamma-direction, 

. 2/iu  2 

in  /lb-sec 

Same  as  PIMA,  only  in  the  beta-direction 

Calculated  pressure  on  panel  (uniform  dis- 
tributions only),  psi 

Pressure  P^,  psi  (Load  option  3) 

Pressure  P^,  psi  (Load  Option  3) 

Output  frequency  - integration  steps  per 
printout 


PRES(3*MGMB) 


Working  array  in  RELAXP,  dimensionless 


*PRT (NTIME , NPY , NPX) 
PRTT(NPY.NPX) 

*PT(20) 

PX(3*MGMB) 

Q1 

Q2 

RA(NBAR.MBAR) 

RFR 

RHO 

*RHOM(NL) 

RR(4) 

RRES (3*MGMB) 
RTRIAL(5) 

SAC (NL) 

SAT(NL) 

*SIGO 

SIG02 

SIGTTl(LMAX) 

SIGX(3*MGMB) 
SIGXT1 (LMAX) 


Pressure  specified  on  panel  vs  time,  psi 
(Load  Option  2) 

Temporary  storage  for  pressure  on  panel 
after  interpolation  on  time,  psi  (Load 
Option  2) 

Table  of  uniform  pressures  specified  on  panel, 
psi  (Load  Option  4) 

Working  array  in  RELAXP,  dimensionless 

Constant  A (Load  Option  1) 

Constant  B (Load  Option  1) 


Range  from  detonation  point,  R,  in. 


Not  used 


2 4 

Density  of  panel,  lb-sec  /in 

2 4 

Density  of  layer,  lb-sec  /in 
Reaction  force  at  corners  of  panel 
Working  array  in  RELAXP,  dimensionless 
Not  used 


Compressive  yield  stresses  for  metal  material, 
compressive  ultimate  stress  for  plastic  material, 
psi;  not  used 

Tensile  yield  stress  for  metal  materials, 
tensile  ultimate  stress  for  plastic  material, 
psi ; not  used 

Yield  stress  for  elastic-plastic  option,  psi 

2 a 

Constant,  equal  to  SIGO  squared,  lb  /in  . 

Stress  component  at  time  of  last  yielding, 

°60  ’ Psi 

Working  array  in  RELAXP,  dimensionless 
Stress  component  at  time  of  last  yielding, 
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SIGXX1 (LMAX) 

SINB(MB*NBAR) 
SING(MG*MBAR) 
SIN2B(MB*NBAR) 
SIN2G (MG*MBAR) 
SIA(NGNBT) 

S2A(NGNBT) 

S3A(NGNBT) 

SAA(NGNBT) 

S5A(NGNBT) 

S6A(NGNBT) 

SMAX 

STT(LMAX) 

SXT(LMAX) 

SXX(LMAX) 

*THETAO 

*THNU (NL) 

*TITLE(20) 

TMAX 

*TNU 

*TPRIME 

*TSTOP 

*TT(20) 


Stress  component  at  time  of  last  yielding, 
°xx’  Psl 

Sines  of  beta  functions,  sin((n+l)£j) 

Sines  of  gamma  functions,  sin  ((aH-1)-^) 
Sines  of  beta  functions,  sin((n)f3j) 

Sines  of  gamma  functions,  sin((m)y^) 

Stress  component,  °™x»  psi 


Stress  component,  a™,  psi 
00 


Stress  component,  am  , psi 

Xu 


Stress  component,  aD  , psi 
xx 


Stress  component,  , psi 
00 


Stress  component,  , psi 
Not  used 


Stress  component,  o , psi 
00 

Stress  component , a , psi 

Stress  component,  o , psi 

Total  angle  subtended  by  cylindrical  panel 
0q>  deg>  or  width  of  flat  panel,  b,  in. 

Poisson’s  ratio  in  the  theta-direction,  v 
dimensionless  ® 

Descriptive  title  of  case 

Not  used 

Poisson's  ratio  for  elastic-plastic  option 
Time  t ' (Load  Option  3) , sec 
Integration  stop  time,  sec 
Time  table  (Load  Option  4),  sec 
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TTNU(LMAX) 


UB(NGNBT) 
UG (NGNBT) 
UU(MB,MG) 
U1 (MB, MG) 


Value  of  v , dimensionless 
s 

Value  of  U,  dimensionless 
Value  of  U , dimensionless 

p 

Value  of  U , dimensionless 

y 

Displacement  coefficient,  U , dimensionless 

tnn 

Displacement  coefficients,  U , representing 

mn 

the  static  conditions,  dimensionless 


V (NGNBT) 


Value  of  V,  dimensionless 


VB (NGNBT) 
VG (NGNBT) 


Value  of  V , dimensionless 
p 

Value  of  V , dimensionless 


VRX(2*NBAR-2) 


VRT (2*MBAR-2) 


Normal  reactive  force  per  unit  length  along 
boundary,  V^,  lb /in. 

Normal  reaction  force  per  unit  length  along 
boundary,  Vx>  lb /in. 

4 

Shock  velocity,  equal  to  5.88x10  in/sec 
(Load  Option  1) 


W(MB,MG) 

VX0(3*MGMB) 


Displacement  coefficients,  V , dimensionless 

mn 

Initial  velocity  coefficients 


VI (MB, MG) 


Displacement  coefficients,  V , representing 

mn 

the  static  conditions,  dimensionless 


W (NGNBT) 


Value  of  W,  dimensionless 


WB (NGNBT) 

WBB (NGNBT) 
WBBB(2*MBAR-2) 
WG (NGNBT) 


Value  of  W , dimensionless 

P 

Value  of  W , dimensionless 
pp 

Value  of  W00 . , dimensionless 
ppp 

Value  of  W , dimensionless 


WGB (NGNBT) 
WGBB(2*NBAR-2) 
WGG (NGNBT) 
WGGB (2*MBAR-2) 
WGGG (2*NBAR-2) 


Value  of  W .,  dimensionless 
YP 

Value  of  W OQ,  dimensionless 
Ypp 

Value  of  W , dimensionless 
YY 

Value  of  W „,  dimensionless 
YY6 

Value  of  W , dimensionless 
YYY 
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WW(MB,MG) 
W1 (MB, MG) 

XB(NBAR) 

XG(MBAR) 

XJ 

XJ2 

XJ3 

XJ4 

XJ5 

XKTT 

XKXT 


Displacement  coefficients,  W , dimensionless 

mn 

Displacement  coefficients,  W , representing 
the  static  conditions,  dimensionless 

Integration  positions  in  the  beta-direction 
inches  for  flat  panel,  degrees  for  curved 
panel 

Integration  positions  in  the  gamma-direction, 
inches 


Constant,  J,  equal  to  180/0^  (dimensionless) 

for  curved  panel  and  ir/b  (inches)  for  flat 
panel 

2 

Constant,  J 
Constant,  J/L 
Constant , 2J 
Constant,  2J/L 

Temporary  value  of  strain,  , in/in. 

uw  • 

Temporary  value  of  strain,  in/in. 


XKKK 

XL 

*XLP 

XLP1 

XLP2 

XL1 

XL2 

XL3 

XU 

XL5 

XL7 


Temporary  value  of  strain,  e^,  in/in. 

Constant,  X./ir,  for  flat  panel  (inches), 
l/va  for  curved  panel  (dimensionless) 

Length  of  panel,  i,  inches 
2 

Constant,  2L  R 

Constant,  2LR 

Constant,  1/L 
2 

Constant,  L 

Constant,  2/L 
2 

Constant,  2L 

2 

Constant,  2L  R 
2 

Constant,  1/L 
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*XP(NPX) 


X-position  of  pressure-mesh  (Load  Option  2) 
in. 


1 


i 


XRES (3*MGMB) 
XX(3*MBMG) 
*XXNU (NL) 

XXI (3*MGMB) 
XIA(NGNBT) 

X2A(NGNBT) 

X3A(NGNBT) 

X4A(NGNBT) 

X5A (NGNBT) 

X6A(NGNBT) 
*YP (NPX) 

YY (3*MGMB) 

ZA(2) 

ZB  (2) 
ZC(2*NL) 

*ZEE 

ZF(LBAR) 

ZG(LBAR) 


Working  array  in  RELAXP,  dimensionless 

Array  composed  of  UU,  W,  and  WW,  dimensionless 

Poisson's  ratio  in  the  x-direction,  v , 

X 

dimensionless 

Working  array  in  RELAXP,  dimensionless 
Strain  component,  e™  , in/in. 

Strain  component,  e™  , in/ in. 

0 0 


Strain  component,  emQ,  in/in. 

X0 


Strain  component,  e“  , in/in. 


Strain  component,  e°,  in/in. 

00 


Strain  component,  e°  , in/in. 

Y-position  of  pressure-mesh  (Load  Option  2), 
in.  or  deg 

Array  composed  of  acceleration  coefficients, 

U , V , W , 1/sec2 
mn  mn  mn 

ZB  normalized  with  a,  ZB/a,  inches  for  flat 
plate,  dimensionless  for  a curved  panel 

Not  used 

Not  used 


Distance  from  panel  to  detonation,  Z,  in. 
(Load  Option  1) 

ZH  normalized  with  a,  ZH/a,  inches  for  a 
flat  plate,  dimensionless  for  a curved  panel 

Gaussian  station  squared,  dimensionless 


I 

1 


. 


] 

1 

■ 
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ZH(LBAR) 


z coordinates  corresponding  to  the  Gaussian 
integration  points  through  the  thickness  for 
the  elastic-plastic  option,  in. 

ZZ1(9)  Not  used 

3.3  Program  Operation 

The  DEPROP  program  is  written  in  FORTRAN  IV  and  consists  of  19  user 
supplied  routines  on  approximately  2700  cards.  The  code  was  developed 
on  the  Control  Data  Corporation  (CDC)  6600  computer  under  the  SCOPE 
3. A. 3 operating  system. 

In  order  to  minimize  the  amount  of  central  memory  core  required  to 
execute  the  program,  the  user  should  eliminate  at  least  one  of  two 
options  prior  to  loading.  This  choice  of  options  is  either  the  elastic 
static  solution  capability,  where  subroutines  RELAXP  and  SOLVE  are 
required,  or  the  elastic-plastic  option  where  subroutines  LIST2  and 
SIGMA  are  required.  These  options  are  outlined  in  Tables  5 and  6, 
and  the  corresponding  program  core  requirements  are  given.  Elimination 
of  the  unnecessary  subroutines  can  be  accomplished  by  removing  them  from 
the  deck  (or  file)  completely,  replacing  them  with  dummy  routines,  or 
using  an  SLOAD  instruction  to  selectively  load  the  required  routines. 

The  use  of  blank  common  enables  the  program  to  load  and  execute  at  the 
same  core  level.  Input  and  output  are  equated  with  logical  files  TAPE5 
and  TAPE6,  respectively,  and  there  are  no  additional  tape  or  disk  file 
requirements  for  operating  the  code. 

The  FTN  compiler  has  been  used  to  compile  the  code  under  "0PT=1" 
and  "R=2"  options.  Compilation  requires  approximately  50000.,  cells  and 

O 

20  seconds  CP  time. 

Computation  times  will  vary  considerably,  depending  on  the  panel, 

the  complexity  of  the  model,  whether  or  not  the  solution  goes  inelastic. 

and  the  response  time  requested.  A very  rough  approximation  for  an 

-4 

inelastic  response  is  5 x 10  CP  seconds  per  mode,  per  Integration  grid 
point,  per  time  step  of  response. 
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HP  P 


- 


Program  Response 
Option 


Elastic  (Dynamic) 


Input  Parameters 


NDERV 


Elastic  (Static)  1 or  3 


Elastic-Plastic 

(Dynamic) 


Static,  followed 
by  Elastic- 
Plastic  (Dynamic) 


*Code  compiled  under  FTN  option  1 


Subroutines 

Eliminated 

Core  Required 
to  Load  and 
Execute* 

RELAXP,  SOLVE, 

125  K8 

LIST2,  SIGMA 

LIST2 , SIGMA 

200  Kg 

RELAXP,  SOLVE 

2!0K8 

None 

262  K8 

TABLE  6.  PANEL  RESPONSE  OPTIONS 


3.4  Program  Input  Data 


i 


Specific  instructions  are  provided  in  this  section  to  enable  the 
user  to  provide  the  necessary  card  input  for  the  DEPROP  program.  Pre- 
ceding these  instructions  are  several  paragraphs  containing  general 
remarks  and  amplification  of  some  of  the  specific  instructions. 

The  input  data  are  specified  in  groups,  where  each  group  begins  on 
a separate  card.  More  than  one  card  may  be  required  for  a group,  how- 
ever. The  variable  type  and  format  corresponding  to  each  data  group  is 
given  in  parentheses  in  the  input  instruction  and  is  always  in  fields 
of  12.  For  convenience,  floating  point  numbers  can  be  left  justified 
the  field  as  long  as  the  exponent  is  right  justified.  Also,  zero  values 
can  be  replaced  by  a blank  field.  Columns  73  through  80  are  not  used 
for  data  and  can  be  used  for  card  identification  or  other  purposes. 

All  input  parameters,  where  appropriate,  should  be  compared  with 
the  maximum  dimensions  provided  for  in  the  program,  as  delineated  in 
Table  3.  This  is  very  important  since  the  program  does  not  attempt  to 
check  the  input  for  such  violations. 

Group  1 provides  for  the  execution  of  several  jobs  (or  cases)  in 
the  same  run.  All  subsequent  data  groups  (2  to  23)  must  be  repeated  for 
each  case. 

The  panel  type,  response  option,  and  debug  options  are  specified  in 
Group  3.  It  is  important  to  check  Table  5 to  insure  the  correct  sub- 
routines are  included  to  the  response  option  selected.  It  is  suggested 
that  the  first  debug  option  (NDBUG=0)  be  selected. 

Group  4 contains  the  number  of  modes  to  be  used  in  the  solution 
and  the  number  of  integration  points  to  be  used.  The  accuracy  of  the 
solution  is  based  on  the  degree  of  convergence  of  stress  and  strain 
quantities.  These  quantities  converge  less  rapidly  than  the  radial 
displacement.  Also,  cases  involving  a clamped  edge  condition  will 
converge  less  rapidly  than  simply-supported  cases.  Since  both  computer 
time  and  accuracy  increase  with  more  modes  and  points,  a trade-off 
usually  becomes  necessary.  It  has  been  found  that  about  25  modes 
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give  an  acceptable  general  solution  for  most  panels,  but  more  modes  are 
required  for  better  accuracy  in  determining  edge  strain  and  reaction 
force  quantities  for  clamped  panels. 


The  actual  mode  numbers  are  specified  in  Groups  5 and  6.  The 
maximum  value  that  the  mode  numbers  can  assume  in  the  program  is  19. 

When  symmetry  is  taken  in  either  direction  (Group  7),  or  if  the  pressure 
loading  is  symmetrically  oriented,  only  the  odd  numbered  modes  (1,  3, 

5,  ...)  are  required  in  that  direction. 

Spatially,  the  desired  number  of  integration  points  (MBAR  and  NBAR) 
for  a full  panel  should  be  approximately  two  times  the  maximum  mode 
number  used  in  that  direction,  plus  three.  However,  when  NBN  or  MGM  is 
large,  this  condition  may  not  be  satisfied  for  nonsymmetrical  panels, 
since  MBAR  and  NBAR  are  dimensioned  at  23  in  the  program  (see  Table  3). 
For  symmetric  solutions,  MBAR  (or  NBAR)  need  only  be  approximately 
one-half  the  value  for  a full  panel  since  only  one-half  (or  one 
quarter)  of  the  panel  is  actually  analyzed  in  the  solution.  For  a 
nonsymmetric  condition,  MBAR  (or  NBAR)  must  be  an  odd  number.  For 
an  elastic-plastic  solution,  a minimum  of  four  integration  points 
through  the  thickness  is  recommended,  and  a maximum  of  six  is  provided 
in  the  program. 

In  Group  8,  the  user  is  given  the  option  of  a purely  elastic 
solution,  or  an  elastic-plastic  solution.  The  elastic-plastic  option 
will  tend  to  be  slower  and  require  more  computer  memory.  It  should 
be  noted  that  honeycomb  panels  are  reduced  to  an  equivalent  single- 
layered panel  for  elastic-plastic  response.  Again,  Tables  5 and  6 
should  be  consulted  to  insure  the  program  is  compatible  with  the  re- 
sponse option  selected. 

Groups  9 and  10  provide  a mechanism  for  selecting  a maximum  of 
49  modal  combinations  from  a 13  by  13  combination  array  (MG“MB*13) . 

Thus,  the  more  significant  modal  combinations  for  an  optimal  solution 
with  respect  to  accuracy  and  computer  time  can  be  selected  and  the 
other  combinations  eliminated.  A general  rule  of  thumb  is  to  eliminate 
the  higher  frequency  modes  which  are  usually  associated  with  modal 
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combinations  having  the  larger  MG+MB  values.  An  example  of  this  would 
be  the  selection  of  MG=MB=7  for  a symmetric  problem,  but  eliminating  24 
combinations  as  indicated  in  Table  7. 


TABLE  7.  EXAMPLE  OF  MODAL  SELECTION 


1 3 5 7 9 11  13 

Gamma  Mode 


The  relative  importance  of  each  modal  combination  can  be  evaluated  by 
examining  the  response  output  and  comparing  the  magnitudes  of  the 
displacement  coefficients. 

Groups  11  and  12  are  responsible  for  selecting  the  points  in  the 
integration  grid  for  which  printout  of  strains,  stresses,  displacements, 
reactive  forces  (or  boundaries),  and  pressures  is  required.  Strains  and 
stresses  are  computed  at  the  inner  and  outer  surfaces  of  the  panel. 

Each  point  in  the  grid  is  designated  by  a pair  of  integers,  the  first 
integer  referring  to  the  gamma-position,  the  second  to  the  beta-position. 

• Actual  positions  are  found  from 

x - j ilrii  x * 1 m 

^ (M— 1 ) 

(symmetric  in  x-direction) 

x **  l I = 1, . . . ,M 

M * (full  in  x-direction) 
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and  similar  expressions  for  y (or  0).  For  example,  the  corner  point  in 
a symmetric  panel  would  be  numbered  (1,1);  the  center,  (MBAR,NBAR) . 

Group  20  contains  the  modal  components,  6 , for  the  Initial  radial 

mn 

imperfections.  The  analyst  must  compute  the  6 ' s from  measured  data 

mn 

using  the  integration  technique  applied  to  Fourier  series  coefficients. 
Generally,  such  data  will  not  be  available,  and  zero  values  should  be 
specified  for  the  6^ 1 s . The  capability  of  considering  initial  imper- 
fections also  enables  the  analyst  to  determine  the  sensitivity  of  panel 
response  to  initial  imperfections. 

Group  21  provides  the  integration  time  increment,  the  response  stop 
time,  and  the  printout  interval.  If  the  user  specifies  a zero  time 
increment,  the  program  computes  an  appropriate  At  which,  in  most  cases, 
will  give  a stable  solution.  Because  it  is  approximate,  the  analyst  may 
want  to  make  comparable  runs  using  different  At's.  In  general,  an 
elastic  solution  which  is  numerically  stable  will  be  accurate.  Hence, 
the  optimum  At  is  the  largest  which  remains  stable.  For  an  elastic- 
plastic  solution,  however,  the  accuracy  of  the  solution  may  deteriorate 
slightly  as  the  point  at  which  the  solution  diverges  is  approached. 

Once  a time  increment  is  selected,  it  should  be  valid  for  moderate 
changes  in  response  level. 

Although  the  stop  time  can  vary  a great  deal,  the  total  number  of 
integration  steps  required  to  capture  peak  response  will  be  roughly 
between  500  to  1500.  One  exception  to  this  may  be  a curved  panel 
experiencing  snap-through  buckling,  in  which  case  considerably  larger 
response  times  may  be  required.  A printout  frequency  of  once  every  20 
steps  is  usually  adequate  for  monitoring  the  response  time  history. 

Groups  22  to  33  provide  for  the  appropriate  pressure  load  on  the 
panel.  The  user  should  refer  to  Figure  6 or  7 for  a definition  of 
certain  input  parameters. 


Group  1:  (112)  NCASES 


Number  of  cases  to  be  run.  (NCASES) 


Group  2:  (20A4)  TITLE 

Identifying  title  per  run,  date,  etc.  Free  field. 
(TITLE) 


Group  3:  (3112)  KTYPE,  KDS,  NDBUG 

Code  designating  panel  type  (KTYPE) 

1,  Single-layer  panel 
3,  Honeycomb  panel  (3  layers) 

5,  Multilayer  panel  (elastic  response  only) 


Response  option  code  (KDS) 

1,  static  only 

2,  dynamic  only 

3,  static,  followed  by  dynamic 

Debug  option  (NDBUG) 

0,  no  debug  output 

1,  most  debug  output 

2,  all  debug  output 


Group  4:  (5112)  MG,  MB,  MBAR,  NBAR,  LBAR 

Number  of  gamma  modes  to  be  used  (MG) 

Number  of  beta  modes  to  be  used  (MB) 

Number  of  gamma  Integration  points  actually  used  over 
the  portion  of  the  panel  analyzed.  Must  be  an  odd 
number  for  full  panel  (see  Group  7).  (MBAR) 

Number  of  beta  integration  points  actually  used  over 
the  portion  of  the  panel  analyzed.  Must  be  an  odd 
number  for  full  panel  (see  Group  7).  (NBAR) 

Number  of  z integration  points  used  through  the  thick- 
ness. (LBAR) 

( Not  needed  for  NDERV=1  (see  Group  8)] 


Group  5:  (6112)  (MGM(I) , 1=1,  MG) 

Gamma  mode  numbers,  m. 
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Group  6:  (6112)  (NBN(l) , 1=1,  MB) 

Beta  mode  numbers,  n. 


Group  7:  (2112)  NSYMG , NSYMB 

Symmetry  code  in  gamma  direction  (NSYMG): 

0,  symmetry  assumed  (0  < y £ ir/2) 

1,  no  symmetry  (0  y £ it) 

Symmetry  code  in  beta  direction  (NSYMB): 

0,  symmetry  assumed  (0  <_  B £ ir/2) 

1,  no  symmetry  (0  B it) 


Group  8:  (3112)  NPLT,  NBND,  NDERV 

Panel  type  (NPLT): 

0,  flat  panel 

1,  cylindrical  panel 

Boundary  condition  code  (NBND) . 


y-direction 

1,  clamped-clamped 

2,  simple-simple; 

3,  clamped-clamped; 

4,  simple-simple; 

5,  clamped-simple; 

6,  clamped-clamped; 

7,  clamped-simple; 

8,  simple-simple; 

9,  clamped-simple; 


B-direction 

clamped-clamped 

simple-simple 

simple-simple 

clamped-clamped 

clamped-clamped 

clamped-simple 

simple-simple 

clamped-simple 

clamped-simple 


Note:  Whenever  a clamped-simple  condition  is  selected,  the  full 

panel  is  analyzed  in  that  direction,  and  NSYMG,  NSYMB,  MBAR 
and  NBAR  should  reflect  this. 


Response  option  (NDERV) : 

1,  elastic  only 

2,  elastic-plastic 


Group  9:  (112)  NNOUT 

Number  of  modal  combinations  to  be  eliminated  from 
solution  (NNOUT) . 

(0  <_  NNOUT  < MG*MB) 


If  NNOUT=0 


skip  to  Group  11. 


Group  10:  (2112)  MOUT(I),  NOUT(I) 

Gamma  mode.  (MOUT(I)) 
Beta  mode.  (NOUT(I)  ) 


Repeat  Group  10  for  1*1, 
In  any  order. 


NNOUT . The  cards  in  Group  10  may  be  arranged 


Group  11:  (112)  NKP 

Number  of  spatial  points  at  which  printout  of  stresses, 
strains,  displacements,  reactive  forces  and  pressures 
are  requested.  If  NKP=0,  all  of  the  above  information 
will  be  suppressed. 

(NKP) 


Group  12:  (2112)  KPG(I),  KPB(I) 

Integration  point  in  gamma-direction  at  which  printout 
is  requested.  Points  are  ordered  1-MBAR,  beginning  at 
Y*0,  and  evenly  spaced  from  there.  (KPG(I)) 

Integration  point  in  beta-direction  at  which  printout 
is  requested.  Points  are  ordered  1-NBAR,  beginning  at 
6*0,  and  evenly  spaced  from  there.  (KPB(I)) 


Note:  These  two  indices  are  taken  as  pairs  where  each  pair  desig- 

nates a particular  spatial  point.  The  pairs  may  be  specified 
in  any  order. 


Group  13:  (112)  NL 

Number  of  layers.  (NL) 

(NL  must  be  1 for  KTYPE-1,  and  3 for  KTYPE=3) 


(THETAO) 


Group  14:  (3F12.1)  XLP,  THETAO,  A 

Full  length  of  panel,  in.  (XLP) 


Full  width  of  flat  panel,  b (short  direction), 
in.  (NPLT=0) 


Full  subtended  angle  of  cylindrical  panel,  0n, 
deg.  (NPLT=1)  0 

Radius  of  cylindrical  panel,  in.  (A) 

(Not  needed  for  NPLT=0) 

If  NDERV=2.  skip  to  Group  18. 


Group  15:  (2F12.1)  HM(I),  RHOM(I) 

Distance  (h)  from  the  inner  panel  surface  to  the  outer 
surface  of  layer  I,  in.  (HM(I)) 

9 / 

Mass  density  of  layer  I,  lb-sec  /in\  (RHOM(I)) 


Group  16:  (5F12.1)  EX(I),  ET(I),  XXNU(I),  THNU(I),  GXT(I) 

Modulus  of  elasticity  in  the  x-direction,  psi.  (EX(I)) 
Modulus  of  elasticity  in  the  theta-direction,  psi.  (ET(I)) 
Poisson's  ratio  in  the  x-direction.  (XXNU(I)) 

Poisson's  ratio  in  the  theta-direction.  (THNU(I)) 

Shear  modulus,  psi.  (GXT(I)) 


Group  17:  (2F12.1)  SAT(I),  SAC(I) 

Tensile  yield  stress  for  metal  panels;  tensile  ultimate 
stress  for  plastic  panels,  psi.  (SAT (I)) 

Absolute  value  of  compressive  yield  stress  for  metal 
panels;  absolute  value  of  compressive  ultimate  stress 
for  plastic  panels,  psi.  (SAC(I)) 

Repeat  Groups  15-17  for  1=1,  NL. 

Skip  to  Group  20. 
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Group  18:  (3F12.1)  HM(I),  RHOM(I),  EM(I) 

Distance  (h)  from  inner  shell  surface 
surface  of  layer  I,  in.  (HM(I)) 

2 4 

Mass  density  of  layer  I,  lb-sec  /in  . 

Modulus  of  elasticity,  psi.  (EM(I)) 

Repeat  Group  18  for  1*1 f NL. 

Group  19:  (4F12.1)  TNU,  SIGO,  EP,  EPSIF 

Poisson's  ratio.  (TNU) 

Yield  stress  for  a metal  panel,  psi.  (SIGO) 

Strain  hardening  modulus  (E^),  psi.  (EP) 

Ultimate  strain,  in/in.  (EPSIF) 

(not  necessary) 

Group  20:  (6F12.1)  ((FG(N,M) , N=1,MB),  M=1,MG) 

Modal  displacement  coefficients  for  initial  radial 
imperfections,  in.  (FG(N,M)) 

Group  21:  (3F12.1)  DELTIM,  TSTOP,  PRINT 

Integration  time  increment,  sec.  If  DELTIM=0.0,  the 
program  determines  the  time  increment  required  for 
stability.  (DELTIM) 

Integration  stop  time,  sec.  (TSTOP) 

Print  frequency  (integration  steps  per  printout).  If 
PRINT=0.0,  printout  of  intermediate  data  will  be  sup- 
pressed. (PRINT) 

If  KDS»2,  skip  Group  22. 

Group  22:  (F12.1)  PS 

Uniform  static  pressure  load,  psi.  Can  be  either 
positive  or  negative  value. 


in  the  outer 
(RHOM(I)) 


If  KDS=1,  skip  Groups  23-33. 


Group  23:  (112)  NLOAD 

Dynamic  load  option 

1,  special  Eglin  analytical  function  over  space 
and  time.  (See  Figure  6.) 

2,  discrete  point  by  point,  time  by  time  distri- 
bution 

3,  spatially  uniform,  with  an  analytical  function 
for  time  history.  (See  Figure  7.) 

4,  spatially  uniform,  with  a discrete  time  history. 


If  NL0AD=2 , skip  to  Group  25 . 
If  NL0AD=3,  skip  to  Group  31. 
If  NL0AD=4 , skip  to  Group  32. 


Group  24:  (2F12.1)  ZEE,  PHI 

Distance  of  detonation  from  panel,  Z,  in.  (ZEE) 

Angle  projectile  trajectory  makes  with  the  normal  to 
the  panel  (z-axis),  ()>,  degrees.  (PHI) 

Skip  Groups  25-33 

Group  25:  (3112)  NPX,  NPY , NTIME 

Number  of  spatial  points  in  the  gamma-direction  at 
which  pressures  are  to  be  specified.  (NPX) 

(Must  be  at  least  2) 

Number  of  spatial  points  in  the  beta-direction  at 
which  pressures  are  to  be  specified.  (NPY) 

(Must  be  at  least  2) 

Number  of  times  specified  in  the  pressure-time  history. 
(NTIME)  (2  < NTIME  < 6) 


Group  26:  (F12.1)  DTIM 

Time  interval  between  samplings  (DTIM).  The  time  history 
for  each  point  has  the  same  time  interval,  but  distinct 
delay  times  (Group  29) 

Note:  Be  sure  to  allow  for  first  point  to  be  engulfed  at  time=0. 

Program  will  extrapolate  data  past  last  time  in  table. 
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Group  27:  (6F12.1)  (XP(I),  1=1,  NPX) 

x-posltions  at  which  time  histories  are  specified, 
in.  (XP) 


Group  28:  (6F12.1)  (YP(I),  1=1,  NPY ) 

y-positions  at  which  time  histories  are  specified, 
in.  (YP) 

Group  29:  (6F12.1)  (DET(J.I),  J=l,  NPY) 

Delay  time  for  pressure  wave  to  reach  grid  point, 
sec  (DET)  (One  point  must  have  delay  time  of  zero) 

Repeat  Group  29  for  1=1,  NPX. 

Group  30:  (6F12.1)  (PRT(K,J,I),  K=l,  NTIME) 

Pressure  for  each  time  and  grid  point,  psi 
(PRT) . 

Repeat  Group  30  for  J=l,  NPY. 

Repeat  Group  30  again,  for  1=1,  NPX. 

Skip  Groups  31-33 

Group  31:  (6112)  PP1,  PPO,  TTO,  TPRIME , AA,  ANN 


Pressure,  p1 , psi  (PP1) 

Pressure,  p^,  psi  (PPO) 

Time,  tQ,  sec  (TTO) 

Time,  t',  sec  (TPRIME) 


Parameter  a,  dimensionless  (AA) 

Parameter  n,  dimensionless  (ANN) 

Note:  See  Figure  7 for  definitions. 

Skip  Groups  32  and  33 
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Group  32:  (112)  NTIME 


Number  of  points  to  be  specified  in  point-by-point 
load  description  (NTIME).  (2  £ NTIME  £ 20). 

Note:  Be  sure  to  include  time=0  and  also  an  end  time  which  exceeds 

TSTOP.  Otherwise,  the  last  value  in  the  table  will  be  used. 

Group  33:  (2F12.1)  TT(I),  PT(I) 

Time,  sec.  (TT(I)) 

Pressure,  psi.  (PT(I)) 

Note:  One  time  and  one  pressure  per  card. 

Repeat  Group  33  for  1=1,  NTIME 

Repeat  Groups  2 to  33  for  each  additional  case,  as  specified  in  Group  1 
3.5  Program  Output 

The  output  for  DEPROP  is  directed  to  the  printer.  Although  program 
output  is  largely  self-explanatory,  the  normal  output  is  described  in 
detail  in  Table  8.  Where  possible,  the  corresponding  program  variable 
is  given  parenthetically. 

Certain  errors,  if  detected  by  the  program  during  execution,  are 
brought  to  the  user's  attention  by  means  of  a printed  error  message. 
Table  9 provides  a list  of  such  messages,  along  with  an  indication  of 
the  routine  associated  with  the  message  and  the  subsequent  action  the 
program  takes.  In  most  cases,  the  program  will  cycle  back  to  attempt 
the  next  case  if  it  cannot  continue  with  the  current  one. 
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TABLE  8.  DEPROP  STRUCTURAL  RESPONSE  OUTPUT 


Time-History  Output 

Time  from  shock  arrival,  sec  (TIME) 

Normalized  axial,  tangential,  and  radial  displacement  modal 
coefficients  for  all  modes,  with  the  beta  mode  index  varying 
most  rapidly  ((UU(J,I),  VV(J,I),  WW(J,I),  J=1,MB),  1=1, MG) 

Table  of  stress-strain  information  for  inner  and  outer  surfaces 
at  each  grid  point  selected: 

X coordinate,  in.  (XG) 

Beta  position,  in.  or  deg  (XB) 

Depthwise  position,  in. 

Axial  strain,  dimensionless 
Circumferential  strain,  dimensionless 
Shearing  strain,  dimensionless 
Axial  stress,  lb/in^ 

Circumferential  stress,  lb/in 
Shearing  stress,  lb/in^ 

Flag  ("*")  indicating  equivalent  strain  has  exceeded  yield 
strain  (elastic  runs  only) 

Counter  indicating  number  of  unloading  and  reyieldings  (KY) 
(elastic-plastic  runs  only) 

Table  of  reactive  force  information  for  each  grid  point  selected 

Normal  reactive  force  (V  or  V ),  lb/in. 

(VRX  or  VRT)  X 9 

Tangential  reactive  force  (N  or  N ),  lb/in. 

(ENX  or  ENT)  X 9 

Reactive  forces  at  corners: 

Reactive  force  (R) , lb 

(omitted  for  panels  clamped  on  all  edges  since  forces  are 
all  zero). 

Table  of  displacement-pressure  information  at  each  grid  point 
selected : 


X-coordinate,  in  (XG) 
Beta-position,  in  or  deg  (XB) 
Axial  displacement,  in  (UF) 
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TABLE  8.  (Concluded) 


Tangential  displacement,  in  (VF) 

Radial  displacement,  in  (WF) 

Pressure,  psi  (PPP) 

Summary  Output 

Message  indicating  whether  run  was  terminated  normally  or  abnormally, 
and  the  time  at  which  computations  stopped,  sec  (TIME) 

Net  CP  time  for  response,  sec  (CPT) 


Number  of  integration  points  which  yielded,  if  any 
(Elastic-plastic  response  only) 


TABLE  9.  ERROR  MESSAGES 


I 


CANNOT  TOTALLY  CORRECT  FOR  OVERSHOOT.  XXX 

An  iterative  process  to  correct  for  overshoot  associated  with 
yielding  has  not  converged  in  five  trials.  This  probably  means 
a numerical  instability  is  creeping  into  the  solution.  Program 
continues  until  such  errors  occur  100  times.  (SIGMA) 

DEPROP  IS  ABORTED  AT  T,  SEC  = XXX. 

DEPROP  cycles  back  to  attempt  next  case.  This  case  is  aborted. 
(DEPROP) 

EPP  IS  OUT  OF  RANGE 

Numerical  instability  detected.  A smaller  At  may  be  required. 
This  case  is  aborted.  (SIGMA) 

IMMEDIATE  RELOADING,  XXX. 

Probable  numerical  instability  creeping  into  solution.  Program 
continues  until  such  errors  occur  100  times.  (SIGMA) 

SINGULAR  MATRIX  IN  S/R  SOLVE. 

The  relaxation  process  has  generated  a singular  matrix  in  deter- 
mining static  equilibrium.  Program  aborts  this  case.  (SOLVE) 

SOLUTION  DIVERGING  IN  DEPROP. 

Very  large  accelerations  have  been  computed  in  DERV2 , indicating 
a numerical  instability.  A smaller  At  may  be  required.  This 
case  is  aborted.  (DERV2) 

SOLUTION  DIVERGING  IN  REI.AXP. 

The  iterative  process  to  find  static  equilibrium  has  failed. 
Program  aborts  this  case.  (RELAXP) 

SOLUTION  IS  UNSTABLE. 

Numerical  instability  has  been  detected  in  elastic-plastic  solu- 
tion. A smaller  At  is  required.  This  case  is  aborted.  (SIGMA) 

THE  VALUE  OF  LBAR  IS  INVALID.  LBAR  = XXX. 

An  incorrect  value  of  LBAR  has  been  specified.  This  case  is 
aborted.  (LEGEND) 


L 
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TABLE  9.  (Concluded) 


TOO  MANY  TRIALS  IN  STATIC  SOLUTION.  MTR  = XXX. 

To  avoid  looping  indefinitely  in  attempting  a solution  representing 
static  equilibrium,  an  upper  limit  of  10  is  placed  on  the  number  of 
trials.  Program  may  need  more  trials  and  adjustment  of  the 
variable  CON  in  RELAX? . Program  cycles  to  next  case.  (DEPROP) 

**  WARNING  **  INCONSISTENCY  IN  SYMMETRY 

A clamped-simple  boundary  condition  has  been  specified,  while  a 
symmetric  solution  has  been  indicated.  Program  continues.  (DSET3) 

**  WARNING  **  - TIME  EXCEEDS  TABLE 

Either  TSTOP  should  be  reduced  as  the  time-pressure  table  extended. 
Program  uses  last  value  in  table  and  continues.  (Load  Option  4) 
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appendix  a 

LIST  OF  ALL  PROGRAM  VARIABLES  AND  THE 
ROUTINES  IN  WHICH  THEY  APPEAR 
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•COMOECK  FIRST 

COMMON  /FIRST/  ICOUNT 
•COMOECK  CNOVA 

COMMON/CNOVA/  CR1U5)  ,DkLUM,6AMMA(41)  , ICUMP,  lNOUT,KALT,KB, 

1 KDAM, KDS,K|RR,KOK, K TYPE , NCALL,*C AS! , NCHPT , N0BU6, NMAS8, NTRI AL 

2 P8 (40), POAM, PPP, PR  I NT, RFR,R TRIAL (5), TIME, TITLE (20)  , TSTOP, 

3 ZZl(9) 

•COMOECK  CLOAO 

COMMON  /CLOAO/  PPl,PPO, TTO,TPRIME,AA, ANN, QT T 1 , 0 T TO, AZ, 

1 JL#  NT IM(,MLOAO, PT (20),TT(20),  ZEE , PMI , U1 , 08, V S, 

2 0ET(10, 10),NPX,NPY,DTIM,PrfT(b, 10, 10) ,XP(10) , TP(10) , 

3 1X1(23), JYJ (23), JLT< 10,10), P«TT( 10, 10), OX  1 (23), 0T1 (23) 
•COMOECK  CSLK1 

C0MM0N/CSLK1/  A,KZ,LBAP,LMAX,MB,MBAR,M6,M0M(13) ,M&MS,MGMB2, 

1 N§AR,NBN(l3),NBN0,NBT,NDEfiV, N6NB , NOT, NPL T,NSYMB,NSYMG»PI 

2 , MUSE ( 1 3, 13) , NB , NS , N6NBT 
•COMOECK  CBIK2 

C0MM0N/CBLK2/  BETR (23) , CC 1 ( 1 3) , CC2 ( 1 3) , CCS ( 1 3) , CCS ( 1 3) , 


CK(b),C0SB(299),C0»6(299),C082B(299) , CQS26 (299) , DPM T , DPR T 1 , 
FPl(299),FP2(294),PP3(299),FP4(i3,2),FP5(299),FPb(299)  , 
FP7(299),FPB(|J,2), 

6AM (23), KC, PIMA (23), PINA (23), SINS (299) , SI N6 (299), 

SIN2B(299),SIN2e(299),XJ,XJ2,XJ3,XJ4,XJS,XL,XLP,XuPl,XLP2, 

RL4>>L2,XL3,XL«,XL5,XL7 

CBlKI 

6X(b),MGQ(S) 


NY2,VX0(147),XX(147), YYC147) 


1 
2 
3 

3 

4 
3 

•COMOECK 

COMMON/ CBLK 3/ 

•COMOECK  CBLK4 

C0MM0M/CBLK4/ 

•COMOECK  CBLK9 

COMMON/CBLK*/  EM(8),ERR(147),F6(13,13),HM(0),MQUT(lb9), 

1 MOOT ( 1 b9) ,RHOM(B) ,01 (1 3, 1 3) , VI ( 13,  13)  , Ml  (13,  13) 

•COMOECK  CBLKb 

COMMON  /CBLKb/  ALT  T (1605), 

AL«t ( 1 BOB), ALXX( 1605), BE1 11  SOS), BE2(1 805), BE3(1»05) , 
EPB(lB0S),ETTt(lS05), 

EXT  1 (1B03),EXX1 (1S05),S16TT1  (1005) ,S16XT1 (1S05),9I6XX1 (1S05 
TTNUOBOB) 

•COMOECK  CBLK7 

COMMON/ CBLK 7/  CNl,CN12,CN13,CN2,CN3,CN4,CN5,CNb,CN7,EL,EP, 

1 EP0,EPP,M,IF1RST, JF IRST , LC, LCMAX, NELP , SI  60, SI602, T NO, TNUSU 


1 

2 

3 

4 
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9U,P(36I),KA(23,23) 


• CUM'JfcC*  CBt** 

C0MM0V/CBLA8/ 

•CUMDtC*  CBIH9 

CUMM04/CBLX9/  BT L (8) , BXL ( 8) , CCRI T (B) , C I M8T ( 3) , ET ( 6) , 

1 EX  18), 8X 1 (B),NLZ(16), NRES, NTIC0,NZP#9AC(8),8AT(B), 8MAX, 

2 TC9IT  (8),  IHNU(8),  TMAX,XXNU(8) , ZC(16) 

•CUMDECX  CBtKIO 

COMMOM/CBUUO/  D*tt(S61),l)*ii(361), 0*0(361),  U(361), 

1 U8(3bl),U6(361),  V(361),VB(361),V6(S6l)»*(361),*B(3bl)  , 

2 488 ( 361 ), N6 (361 ),N6B( 361) **66(361) 
tCUMDECK  CBUK11 

COMMOM/CBUU/  CMU,CH12,CM22,CM33,OM11,OM12,OM22,OH33,FM11, 

1 FM12,FM22,FM33 

•CUMUECK  C8LM2 

COMMON)  /CBUK12/  OEl.*  ( 1 47  ) , IP ( 147 ) , PRES  ( 1 47 ) , PX  ( 1 4 7)  , RRE8  ( 1*7 ) , 

1 SI6X(147),XRI8(147,147),XX1 (147) 

ftCOMOECK  CBI.K13 

CUMMOM/CBL* 1 3/  DC#EC,EP8IF,GC,HBAR,NL#NMOUT,RHO,THETAO 
•COMOECA  C8CX14 

COMMOM  /CBLK14/  NBU9E (23 . 23) , NRC * C 1 , C2, C3, C4, C5, C6, C7 , 

1 OELK, CELT,  4666 (44) , *088 (44) , wfiGB (44) , *688(44) , 

2 VRX(44) , VRT(44),RR(4),ENX(44) , ENT (44) , NKP, KPG (46) , KPB (46) 
•COMOECK  CBlANK 

COMMOM  CN1 0, CN1 1,CN8,CN9,EPB0( 1805), ETT, EXT, EXX, 

1 J NZ (2) , XSUMA (361),KY(1805),NUSE(23,23), 

2 9TT (1805) ,9XT (1805) , 

3 3XX(1805),UU(13,13),VV(13,13),Mrt(13,l3),XB(23),XS(23)rXKTT, 

4 XKXT,XKXX,XlA(361)»X2A(361),XiA(361),X4A(361),X5A(S61)r 

5 Xb A ( 36 1 ) , Z A ( 2 ) , ZB(2) , ZF (6) , Z6(6) , ZM  (6) t 81 A (361 ) , S2A (361 ) , 

6 93A(361),84A(361),95A(361),86A(361) 


• DECK  DCPRQP 

PROGRAM  Of PROP  l INPUT , OUTPU I , r APES* 1 NPUT , T AP£6»0UT PUT ) 

•call  first 

•CALL  CNOVA 
•CALL  CLOAO 
C 

1 FORMAT (6112) 

2 FORMA T (6F 1 2* 1 ) 

3 FORMAT  (20A4) 

NCASf  • 0 
INOUT  • 1 
ICOMP  • S 

ICOUNT  • 37777000000000Q00000B 
RCAO(Sfl)  NCA8EB 
100  Rf AD  (3,3)  (TJTLE(I) , 1*1,20) 

NC  ARE  > NC  ARE  ♦ 1 
KERR  ■ 0 
MTRUL  • 0 
ADAM  ■ 2 

READ  (5,1)  KTVP(,KOR, NOBUG 
IFIINOUT.EB.O)  GO  TO  1400 
NRITiUiSOOO)  (UTL«U),I*1,20) 

GO  TO  (300,400,500,600,700),  ATYPE 
300  4RI TE (6, 3500} 

60  TO  1050 
400  ARITEU,  3600) 

GO  TO  1050 
500  MR  I T C (6, 37 00 ) 

GO  TO  1050 
600  NRI  T(<(6,  SBOO ) 

GO  TO  1050 
700  MR1 T£ (6, 3000) 

60  TO  1050 

1050  60  TO  (1100,1200,1300),  ADS 
1100  NRXTI(6,4300) 

GO  TO  1400 
1200  MRITf (6,4400) 

GO  TO  1400 
1300  NRI TE (6, 4500 ) 

1400  NCALLI  • 2 


DEPROP 
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CALL  PROP 

IF (HERR.GT.O)  GO  TO  1600 
NCAUU  * 1 
CALL  PINIT(O) 

CALL  PROP 

IF(K0S.EQ,1)  GO  TO  1600 
IF (KERR«GT ,0)  GO  TO  1600 
NCALL*  * 0 
HOK  * 0 
CALL1  PlNIT(l) 

RTR1AL(1)«1.0 
lbOO  iMTR  1 AL  • NTR1AL  ♦ 1 
CALL  PROP 

1 60 0 IF(NCASE.LT.NCASES)  GO  TO  100 
1700  STOP 
C 

3000  FORMAT  (1H1V30XVUH0  E P R 0 P//1X,20A4) 

3500  FORMAT  ( 2SH 031 N6LE "LAYER  M|!T AL'  PANEL  ) 

3600  FORMAT  (30H03IN6LE-LA YER  PLASTIC  PANEL  ) 

3700  FORMAT  (25HOHONE VCOMS  METAL*  PANEL  ) 

3800  FORMAT  (27HOHONE VCOMS  PLASTIC  PANEL  ) 

3900  FORMAT  (29M0MULTI-LAYER  PLASTIC  PANEL  ) 

4300  FORMAT  121H0STAT1C  SOLUTION  ONLY) 

4400  FORMAT  (22HOOVNAMIC  RESPONSE  ONLY) 

4500  FORMAT  (37H03TAT1C  SOLUTION  AND  DYNAMIC  RESPONSE) 
END 


DEPROP 
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• DECK  SOLT 

SUBROUTINE  BOLT 

this  subroutine  sets  up  « mode  shares  for  boundary  conditions 
SELtCTEO. 

call  CBLR1 
call  C»L*2 

0IHCNSI0N  C01(20),C02(20) , CDS (20) , CD4 (2o) 

OJNCNt ION  COL (20) # C0A(20) 

DATA  CD  1/1,5056167  314,2, <1997  52*7005,3 ,50 00 1067945,  4.49999953647, 

1 5,50000001994,6,49999999915,7.5,8,5,9.5,10.5,11,5,12.5,13.5, 

2 14.5# 15. 5, 16. 5, 17. 5, IB. 5, 16. 5,  20. 5/ 

DATA  CD2/0, 9625022 14568, 1.00077731109,0.999966450124, 

1 1.00000144969,0,999999937336, 1.00000000270,0.999999999661, 

2 13*1,0/ 

DATA  CD3/ 1.24967633505# 2. 24999976925, 3. 24999999959, 4, 25, 5. 25, 6. 25, 

1 7 .15# 6. 25, 9, 25# 10. 25, 11. 25, 12. 25# 13. 25# 14, 25, 15. 25, 16. 25, 

2 17,15# 16.25# 19.25,20.25/ 

DATA  CD4/ 1.00077 731 192# 1.00000144989, 1,00000000269, 17*1.0/ 


FAC  • SORT (2.0) 

00  100  !•! # 4 
100  CMI)  ■ FAC 
II  « 0 

60  TO  (500,700,500, 700,900,500,900# 700,900) , NSNO 

CLAMPED  • CLAMPED,  6 AMMA . 

500  00  520  1*1 , MG 
M • M0M ( I ) 

COL (1)  • C01(M) 

520  CDA(l)  • CD2 (M) 

540  00  600  M«1,MS 
XI  • CDl(M) 

X2  * COA(M) 

DO  600  1 * 1 , NG I 
II  • II  ♦ l 
X3  * X 1 *6AM ( I ) 
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tXl  = ExP(X3) 
tX*  s tXP(-X5) 

$L  = S1NIX3) 

CL  • COS ( X 3 J 

FPHll)  • -CL  ♦ X2*8L  ♦ . 5* ( 1 , *X2) *E  X2  ♦ .5* ( 1 ,-X2) *EX1 
FP2(II)  » X 1 * (SL  ♦ X J*CL  • .S«(1,*X2)*CX2  ♦ , 5* ( l . *X2) *£x 1) 
FPS(II)  « xl**2« (Cl  • X2*3L  ♦ . 5* ( 1 , »X2) *EX2  ♦ . 5* (1 . -X2) *EX 1) 
IF  U.EU.l)  FP4  ( M,l)  ■ -2,0*X2*X1**3 

IF  (l.ee.MSARJ  FP4(  M,2)  z X1«*3»(-6L  • X2*Cl  ♦ ,5*((1.  • X2)* 
1 EX  1 • (1.  ♦ X2) *EX2) ) 
bUO  CONTINUE 

CX  ( 5)  * I, /FAC 
GO  TO  1000 

SIMPlt  - SIMPLY#  GAMMA , 

700  00  600  M«1,MG 
XI  ■ MBM(M) 

00  600  1*1, NGT 
II  » II  ♦ 1 

X2  * X l *6 AM ( I ) 

X 3 * 8IN(X2) 

FPl(II)  • X 3 
FP2 ( I I ) ■ X1*C0S(X2) 

FPJ(II)  • -X1«*2*X3 
IF  (I. EG. I)  FP«(M,1)  a -X 1**3 
IF  (I. EG, MBA#)  FP<*  (Mf  2 ) » - X 1 *• 3*CUS ( Hi) 

BOO  CONTINUE 

CK(5)  * FAC 
GO  TO  1000 

C CLAMPED  • SIMPLY.  GAMMA, 

900  00  920  1*1 # MG 
M * MGM l 1 ) 

COL  Cl)  * COS (M) 

920  CO A ( 1 ) s CD4 (M) 

GO  TO  540 


1000  II  * 0 
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60  TO  (1100,1300, 1100,1100,1100, 1500, 1300, 1500, 1500),  N0ND 

CLAMPED  • CLAMPED,  BETA, 

1100  00  1120  Ul,Mg 
N • M«N(I) 

COLU)  « COt(N) 

1120  COA(I)  ■ CD2 (N) 

1140  00  1200  N«1,M0 

XI  a COL(N) 

X2  ■ COA(N) 

00  1200  J«1 , NOT 
11  ■ 11  ♦ 1 
X 3 a X1*6ETR(J) 

EX  1 8 EXP (M3) 

EX2  « ExP(-XJ) 

•L  • SlN(Xl) 

CL  • C0S(X3) 

PP5UI)  • *CL  ♦ X2*8L  ♦ ,b*(l ,*X2)*EX2  ♦ ,5*  ( 1 ,-X2)  *EX1 
FP6U1)  • XI*  (SL  ♦ X2*CL  • , 5* ( 1 , *X2)  *EX2  ♦ , B*  ( 1 . -*2 ) *E  X l; 

FP7 (111  • X 1**2* (CL  • X2*SL  ♦ ,5«(1.«X2)*EX2  ♦ ,5* ( 1 , -X2 ) •£» 1 ) 
IF  (J.EO.l)  FPS(  N, 1 ) * •2.0*X2*X1**3 

IF  (J.EQ.NSAR)  FP«(  N,  2)  ■ Xl**3*(-$L  • X2*CL  ♦ ,5*((1.  • X2)* 
1 1X1  • (1,  * X2)*EX2)) 

1200  CONTINUE 

CM6)  * 1,/FAC 
60  TO  1600 

•IMPLflT  • SIMPLY,  beta. 

1300  00  1600  Mai, MB 
XI  a MBM(N) 

DO  1400  Jai , MBT 
II  a II  ♦ 1 
X2  a X 1 *BE  T R ( J ) 

X3  a SINCX2) 

FP5(ll)  ■ X 3 

FP6 ( II)  ■ X 1 *COS ( X 2 ) 

FPT(II)  a -X1**2*X3 
IF  (J.EO.l)  FP8  (M,  1 ) a - X 1 * *.  3 


IF  l j , fcu , NBA* ) F P0(N,2)  s • X 1 * * J*CQ8 1 X2) 
140  j CONTINUE 

CK(b)  » FAC 
£0  TO  1800 

C 

C CLAMPED  • SIMPLE,  SETA, 


15uO 

00  1520 

I«l» 

MB 

N • N8N 

U) 

CDL(I) 

• CDS 

(N) 

1520 

coach 

• C04 

(N) 

GO  10  1 

140 

1 b o 0 

METUPN 

END 
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•DECK  OEftVi 

SUBROUTINE  DERV2 
C 

•CALL  CBLK1 
•CALL  C»LK< 

•CALL  CSLKS 
•CALL  cblko 
•CALL  CSLKS 
• CALL  CBLMO 
•CALL  CBLK11 
•CALL  CBLK14 
•CALL  CNOVA 
•CALL  CBLANK 

c 

IF (NCALL.EO.O)  CALL  PRESS 

C 

HI 

DO  100  Mai, MB 
00  100  Nal,MB 

XF(MU8E(N,M).E8.0)  SO  TO  100 
UU(N*  M)aRX (1) 

VV(N,M)aXX(M6MB*I) 

NN(N«  N)a  XX(M6MB2*I) 
xa|»l 

100  CONTINUE 
Kao 

00  500  Ial,NST 
II  a 0 

IF  (I.E8.1)  II  a i 

IF  (I.EQ.MBAR.AND.NSYMG.EQ.l)  II  a 2 
00  500  Jal , NBT 

IF (NUOE (J,I),EO(0)  BO  TO  500 
NBC  a NBUSf ( J, 1) 

JJ  « 0 

IF  (J.E8.1)  JJ  a l 

IF  (J*E8.NBAR,AND.N8YMB.E0.1)  JJ  * 2 

Kaxtl 

SSlaO.O 

OBtaO.O 

BtSaOfO 


SS4*0,  0 

•85*0.0 

586*0.0 

557*0.0 

558*0,0 

559*0.0 

5810*0.0 

5811*0.0 

5512*0,0 

8813  ■ 0. 

8514  * 0. 

8815  * 0. 

8816  * 0. 

00  400  M*l,M(i 
MM*  (M»1)*NGT>I 
SMaSIMC (MM) 

CM*C086(MM) 

SM2*51N28(MM) 

CM2*C0S26(MM) 

T 1 * FPl(MM) 

T 2 ■ PP2 (MM) 

T3  « PP3(MM) 

T 4 * 0. 

IP  (II.6T.0)  T4  ■ PP4(MfH) 
81*0.0 
53*0.0 
54*0.0 
56*0.0 
57*0.0 
59*0,0 
511*0,0 
514  * 0. 

00  200  N*1,MB 

IP (M08E (M.M).EU.O)  60  TO  200 
MM*  (M«1)*MBT*J 
T5  « 0. 

IP  (JJ.GT.O)  TS  * PPB(N.JJ) 
SM*SIMS(MN) 

CN*COSB (NN) 

SM2«SIM2B(NN) 
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CN2*CQ8281NN) 

UNN«UU(N,M) 

9MN*99(N,M) 

4MN*44 (N*  M) 

91*91+UMN*8N2 
33*8S»CC2(N)*UMN*CN2 
94*84*VMN*8N 
99*99*CC9(N)*VMN*CN 
87  * 4MN*FP5(NN)  ♦ 37 
$9  • 4MN*FP9(NN)  ♦ 39 
911  • 4MN*PP7(NN)  ♦ 311 
91*  • 4NN*T5  ♦ 914 
200  CONTINUE 

9lll91«IN«|9l 
392*81*CC5(M)*CM  ♦ 382 
991*9 J*9M  ♦ 933 
9lllM»9M2t89( 

999*94*CC1 (M)*C*2*S35 
899*89 *9M2+39b 
987  * S7*T  1 ♦ 897 
999  ■ 87*T2  ♦ 888 
989  • 89«T1  ♦ 999 

9910  • 87*T3  ♦ 9910 

9911  • 91 1 *T 1 ♦ 9911 

9912  • 89* T 2 * 3812 

IF  (NiC.EQ, 0,04, NBC . GT  , 1 00 ) SO  TO  400 
8313  « 37*T4  ♦ 3313 
9914  ■ 31 4*  T 1 ♦ 8814 
9815  • 31 1 • T2  ♦ 8915 
9919  • 39*T3  * 8919 
400  CONTINUE 
J(4)*981 
U9(K) *992 
U9  (4 ) *991 
V(K)*99« 

96(419999 
9900*999 
4(41*997 
46(4) *899 
49 (4) *899 
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aGGU)«8S10 
4BB  1*0  BS81 1 
*G8(X)»3912 

IF  (NBC.EQ.O.OR.NBC.GT.IOO)  60  TO  500 
NBC  • I ABB ( NBC ) 

IF  (II.EO.O)  GO  TO  450 
NGCG(NBC)  * 3813 
NGBBtNBC)  ■ 8815 
GO  TO  500 

450  RBBB(NBC)  > $314 
rtGGB(NBC)  « 8814 
50U  CONTINUE 

L 

C CONFUTE  STRAINS  ANO  STRESSES 

C 

K»0 

00  700  I«1,NGT 
00  700  J*1 f NBT 

IF(NJSE(J,I),EQ,0)  GO  TO  700 

K«K*1 

UFtUOO 

UGFbUG(K) 

U8F«UB(K) 

VF»V(K) 

V6F«VG(K) 

VBFsVB(K) 

NFatiU) 

NGFcNS(k) 

nbfrnbck) 

OnGF»ON6(K) 

ONBFaONB(K) 

EXX»XU*(UGF«XU*(NGF*0NGF«0.5*(NGF**2^V6F*»2#UGF**2))) 
ETT«XJ*(VBF*XJ»(WBF*DwBF*Ot5*lN0F**2*VBF**2*UBF**2))) 
EXT>XJ*UBF*XL1*  (VGF*  (1 ,0*XJ*VBF)*KJ*(4GF*(4BF'tOWBF)  ♦D*6F*NBF«UBF* 
1UGF ) ) 

AC  ■ X J*V8F  ♦ Xll*U6F  ♦ 1,0 
«6GF  > N60(K) 

4BBF  • «BB(K) 

M6BF  ■ «GB(K) 

XX X X ■ X17*WGGF*AC 
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XKTT  • XJ2*NBBF*AC 
XKXT  B XJS*WQBF*AC 
IF(N>>tT,lU.O)  60  TO  600 

ETTalf  T-*F*(l,0*XJ*VBF»0.S*i»F)*VF*tKJ*NBF*0.5*VF) 

EXT  »tXT*XLl*(a6F*VF*VGF**F) 

XKTT«XKTT*XJ4*VBF*Xll*UGF.NF 

XKXTaXKXT*XLl*V8F 

AC  • XJ*NBF  ♦ VF 

XKXX  • XKXX  • X17*»I66F*MF 

XKTT  • XKTT  ♦ X J*V8F* (X J*VBF  • *F)  • XJ2»*F**8BF  ♦ 

1 (XJ*VBF  • NF) **2  ♦ AC* ( AC  ♦ X J*NBF ) 

XKXT  a XKXT  • XJ5»*8BF*NF  ♦ XL3*NGF*AC 
BOO  IF  (NOERV.CO.2)  60  TO  640 

S1A(K)  • CMllftCKX  ♦ CM12*ETT  ♦ FMlt*XKXX  ♦ FM12*XKTT 
•<A(X)  • CM22*CT  T ♦ CM) 2*EXX  ♦ FM22*XK T T ♦ FM12*XKXX 
83A(R1  a CM33*CKT  * FM33*XKXT 

B4A  (K)  • DM1 1*KKXX  ♦ DM12*XXTT  ♦ FMU*EXX  ♦ FM12*ET  T 


85A (K)  a DM22*XKT  T ♦ 0M12*XKXX  ♦ FM22*ETT  ♦ FM12*EXX 
S6A (K)  ■ DM33* XK X T ♦ FM33*EXT 

60  TO  660 

640 

CALL  SI6MA  (J,I,K) 

| • 660 

XI A (KJ  • EXX 
X2A  (X)  • ETT 
X3A (K)  • EXT 
XOA(R)  * XKXX 
XSA ( K ) a XKTT 
X6A (K)  • XKXT 

700 

CONTINUE 

IF (KCKR*6T «0)  60  TO  2200 

C 

K « 0 

00  750  1st t N6T 
DO  750  Jai»NBT 

IF  (NUBiWrD.BT.O)  K a K ♦ 1 
IF  (NBUBI(J'I).NE.O)  CAU.  *CIT(X,J,K) 

750 

CONTINUE 

c 

• 

IF(NCALL.EO.l)  60  TO  900 
KZaO 

IF  (K0AM(LT.2,ANDvKC.E0,10)  K2  ■ 1 

■ 
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IF  (PRINT. EQ.O.)  GO  TO  800 
IF(TIME.LT.DPRT)  GO  TO  BOO 
Kl9Kl*l 

IF  (KOAM.LT. 2)  KZ  a 3 
QPRTaDPHT^DPRTl 

C PRINT  RESULTS  ANO/OR  CHECK  MAXIMUMS 
800  IF  (NOERV.EQ.l)  CALL  L1ST1 
IF  (NDERV.E0.2)  CALL  LIST2 
IF  (K2.E0.1.0R.KZ.CQ.S)  KC  a 0 
KC  ■ KC  ♦ 1 
C 

c 

900  I Z»0 

00  1800  IRal,MG 
MMQa(IR»l ) *NGT 
00  1800  I Sal t MB 

IF(MU8E(I8,I«).EU.O)  GO  TO  1800 

NNOa (18*1) *NBT 

1 Z« I Z ♦ 1 

8UR8«0.0 

SVRSaO.O 

8NRS>0 . 0 

KaQ 

00  1700  I«UNGT 
MMSMMO# I 
SMsSING (MM) 

CM1CO8G (MM) 

8M2»SIN2G(MM) 

CM2*C082G (MM) 

Ti  « FPl(MM) 

T 2 * FP2(MM) 

TS  « FPS(MM) 

SU  * 0, 

8V  a 0. 

8p  a 0. 

PRLM  • PIMA(I) 

00  1800  J«1,NBT 

IF  (N08E(J,I).EQ.O)  GO  TO  1800 
K a K ♦ 1 

IF  (NU8E(J,I).EQ.l)  GO  TO  1800 


PftlN  « PINA ( J ) 

NNbNNO* J 
SNstlNS(NN) 

CNsCOSB(NN) 

SN2*iIN2B (NN) 

CN2«COS2B(NN) 

UFbU(X) 

U6FBU0U) 

UBFbUB(K) 

VFbVU) 

VCFbVS(K) 

VBFbVB  (K) 

NFbN(K) 

n6FbNB(K) 

NBFbaBIK) 
rtCSF  ■ N660O 
NMF  s MSB (K ) 

N6BP  B *6»(K) 

OnSFbOmS(K) 

ONBFbONB(K) 

IF(NU.Ifi.O)  PPPBP(K) 

PUb8M*BN2 

PUSbCC5UR)bCH*8N2 
PUBbCCC(IB) *8M*CN2 
PVsBNtABN 

PVGbCC11H»)*CM2*SN 
PVBBCC6a8)*6H2*CN 
PN  B T1*PP5(NN) 

PNC  » T2»FP5(NN) 

PnB  B TIbFPB(NN) 

PW06  a T1*FP5(NN) 

PwBB  B T1*FP7(NN) 

PnSB  B T2*FP6 (NN) 

1100  PEXXUBXLt*PUC*U.0*Xll*U6F) 

PE X XVBXU 7 a VCFbPVC 
PEXXN«X(.7bPN6«(N6F«0N6F) 

PETTI/bX  J2*UBF*PUB 

PETTV  a X J*PVB* (1,0  ♦ XJ*VBF) 

PfcTTN  B X J2*PnB* (NBF  ♦ DNBF) 
PEXTUBXJ«CPUB«U.0*XLI*UGF)*XU1*UBF*PU(») 
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PExfV  • XL  1 * (PVGa  11,0  ♦ XJaVBF)  ♦ XjavCFaPVB) 

PtXTrt  s XJ3a(PM8a(ftCF*DPGF)  ♦ PPGa (aBFaDPBF) ) 
PKXX*aXL7aPPGG 
PKTTpsX J2aPPBB 
PKXTM«XJS*PN6B 
AC  • XJ*VBF  ♦ XL1*UCF 
PKXXU  « XL1*XL7*MG6F*PUG 
PKXXV  • XJ*XL7*MGGF*PVB 
PKXXM  « PKXXP  ♦ XL7*PM6G*AC 
PKTTU  • XJ2*XL1*MBBF*PUG 
PKTTV  • P VBaX JaX J2*PBBF 
PKTTP  a PKTTP  ♦ PmBB*XJ2«AC 
PKXTU  a XJ4aXL7a*G8FaPUG 
PKXTV  * XL3*XJ2*N6BF*PVB 
PKXTP  a PKXTP  ♦ XJS*PW6B*AC 
S1«0.0 
8210.0 

IFCNPwT.EQ.O)  GO  TO  1200 
PETTVaPETTVaPVa(VFaXJapBF)  • XjapFaPVB 
PETTp»PETTp  - PP*(1.0+XJaVBF»PF)  ♦ XJaVFappB 
PEXTV«PEXTV  ♦ XL1*(M6F *PV»PFaPVS) 

PEX  T papf X TP  ♦ XL1*(VF«PmG*V6F*PM) 

PKTTU  * XL1*PUG  ♦ PKTTU 
PKTTV»XJ4aPVB  ♦ PKTTV 
PKTTflaPKTTrt  • PP 
PKXTV«XL1*PVG  ♦ PKXTV 

PKXXP  s PKXXP  • XL7*(P*GG*4F  ♦ PP*P6GF ) 

PKTTV  « PKTTV  ♦ X JaPVBa (4, ax JaVBF  • 3,*PF|)  ♦ 

1 PV*(3.*XJ**BF  ♦ 2.*VF) 

PKTTP  ■ PKTTP  ♦ X J2*PpB8* (»PF  ♦ XL1*UGF)  ♦ Ppa(2,aPF  • 
1 3. *X JtvBF  * XJ2#nBBF)  ♦ PPBax Ja ( 4. ax JapBF  ♦ 3,aVF) 

PKXTV  « PKXTV  ♦ XL3*PV*4GF 

PKXTP  a PKXTn  • XJ5a(PPGB*PF  ♦ PpapGBF  • PPBaPGF)  ♦ 

1 P«G*XL3*(XJ*PBF  ♦ VF) 

81  a DPO(K)  ♦ WF 

82  • VF 

1200  PU  * XLP2*PPP*PU*(«GF  ♦ DpGF ) 

PV  a XLP1 aPPPaPV* ( X Ja (PBF  ♦ 3 PBF)  ♦ 82) 

PP  a x^Pl aPPPaPPa (81  • XLlauSF  • XJaVBF  • 1.0) 

IF  (MOERV.EQ. 1 ) GO  TO  1260 
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1290 


C 

1300 


uoo 


J1  * 16AR*(K*1) 

KtUM  • KSUMA(K) 
IFUSUM.LT.lBARj  CO  TO  1300 


81 

a 

S1A(K) 

62 

• 

824(10 

83 

■ 

S3A (X) 

64 

a 

64AOO 

•9 

■ 

ISA(K) 

60 

• 

St A (X  ) 

Fl 

a 

PCXKUaSl  ♦ PCTTU*62 

♦ 

PEXTU*63 

F* 

a 

F«XXU*64  ♦ PKTTU*BS 

♦ 

PKXTU*6t 

F3 

a 

PCXXV«61  ♦ PCTTV«62 

♦ 

PEXTV*63 

F4 

a 

FKXXV*64  ♦ PKTTV«6S 

♦ 

PKXTV«6b 

F5 

• 

PiXXMtCl  ♦ PETTM*62 

♦ 

PEXTN*6S 

F* 

a 

PKXXa«64  ♦ PKTTa«65 

4 

PXXTN*6b 

FU 

a 

CNlOaFl  ♦ CN1 1*F2 

FV 

CNie*Fl  ♦ CNU«F4 

FN 

* 

CN10*F3  ♦ CN1 l*Ft 

60 

TO  1500 

TOTUMiO.O 
TOTVMaO.O 
TOT  AM  a 0,0 
TOTUB  a 0,0 
TOT  VB  a 0,0 
TOTUB  • 0,0 
DO  1400  KKstfLiAR 
l a J1  ♦ XK 

51  • HBO(KK) 

52  4 SX (KK ) *51 
SI  a SXX(l) 

62  a 8TT(U 

53  a SXTCU 

TOTOM  • TOTUM  ♦ SI •( PCX XU* Cl 
TOTUS  a TOTUB  a 8|*(PKXXU*61 
TOT VM  a TOTVH  * Slft(PEXXV*61 
TOT  VS  a TOTVB  ♦ X2*(PKXXV«61 
TOTaa  a TOTNM  ♦ S1*<PCXXN«61 
T0T4BiT0T*8482*tPKXXWt61  *PK 
FU  a CNSaTOTUM  ♦ CN9*T0TUB 


♦ PETTU*82  ♦ PEXT0*63) 

♦ PKTTU*62  ♦ PKXTU*63) 

♦ PETTV*62  ♦ PCX  TV*63) 
4 PKT T V*62  ♦ PKXTV*63) 

♦ PCTTM*82  ♦ PEXTa»63) 
TTM*62  ♦PKXT«*63) 
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FV  * CN8*TUTVM  ♦ CN8»TUTVB 
F H m CN8*T0TnM  ♦ CN9*T0TMB 


C 

1500  3U  » 3U  ♦ (FU  ♦ PU)*PWIN 
8V  * 3v  ♦ (FV  ♦ PV)*PKIN 
9w  « 8*  ♦ (F*  ♦ P*)*PRIN 
lbOO  CONTINUE 

SURS  ■ SUPS  ♦ PRLM*8U 
8VR3  • 8VR8  ♦ PMLH*8V 
SmRS  • 3RR8  ♦ PRLM*8n 
1700  CONTINUE 

IF  (ABS(SWR8)tST,1.0E30)  60  TO  2150 
YY(IZ)  * -9URS*CKU)*CK(2) 

YY(MCMBMZ)  * •SVR8*CK(|)*CR(«) 
YV(M6MB2«IZ)  * -8*RS«CM5)«CK(b) 

1800  CONTINUE 

60  TO  2200 
C 

2150  KERR  « 1 
ARITE  lb, 2151) 

2151  FORMAT  (30H080LUT ION  DIVCR61N6  IN  OEPROP  ) 
C 

2200  RETURN 
END 
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• DECK 

OOET1 

SUBROUTINE  OBIT! 

C 

• CALL 

CBLK1 

•CALL 

CBLAi 

• CALL 

COLA! 

• CALL 

CBLKA 

• CALL 

C0LR5 

• CALL 

CBLK7 

• CALL 

CBLKB 

• CALL 

CBLA9 

• CALL 

CBLR10 

• CALL 

CBLAI I 

• call 

cblais 

• CALL 

CBLKI4 

•call 

CNOVA 

• CALL 
C 

CBLANK 

C INPUT  OATA 

RfeA0(5*7000)  NS*NO*NBAR,NBAR,LBAR 
READ  (5*7000)  (MSN(l ) , Xa| , M6) 

READ  (5*7000)  (NBN(I ) , 1*1 , MB) 

READ  (5*7000)  NSVMB* N8TMB 
REAO (5* 7000)  NPLT, NBND* NOERV 
REAO (5*  7000)  NNQUT 
IFiNNOUT ,EO,0)  SO  TO  70 
00  SO  1*1  * NNOUT 

SO  REAO(S»  7000)  MOUT(I),NOUT(I) 

70  RCA0(5* 7000)  NKP 

IF  (NAP. tQ, 0)  60  TO  90 
00  00  ls|*  NKP 

60  READ (5* 7000)  KPB(I) ,KPB(I) 

90  IPiKDAM.EB. 1. AND.KT VPE.EQ. 1 ) NDERVaS 
IF  (ROAM.E0. I.AND.KTYPE.EO.I)  NOERV  « 2 
IF  (NOERV, CO. I)  (.BAR  • 1 
REAO  (5*7000)  NL 
IF  (KlYPl.LT, 5)  NL  • 3 
IF  (KTYPK.LT.l)  NL  • 1 
REAO (5*  7 1 00)  XLP*THCTA0*A 
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IFCNPwT.fcU.O)  A»1.0 
IF  (NOERV.EQ.l)  CO  10  ISO 
00  100  1*1, Nl 


100  RiAO(S,7100)  MM(I),RM0M(I),EM(1) 

READ  (5,7100)  TNU,6ICQ,EP,EPBIF 
CO  TO  190 
ISO  00  160  !•!, ML 

READ  (5,7100)  MM(I),RHQM(I) 

REAO  (5,7100)  EX(I),ET(I),XXNU(I),THNU(I),CXT(I) 

160  REAO  (S,7100)  SAT ( I) , 8AC  ( I ) 

IF  (KTYPE,NE.J.AN0,KTYPE.NI,4)  CO  TO  190 
IF  (KOAM.NE.O)  CO  TO  190 
READ  (5,7100)  EC,CC,DC 
190  REAO  (5,7100)  ( (FC ( I , J) , I«1 , MB) , J« 1 , MC) 

REAO (5, 7100)  DELTIM,T8T0P, PRINT 
IF(INOUT'EQ.O)  CO  TO  2100 
C PRINT  OUT  THE  INPUT 
WRITER 6, 7170) 

WRITE  (6,7200)  MC, MB, M8AR, NBAR,LBAR 
WRITE  (6,7210)  (MCM ( I ) , !■ 1 , MC) 

WRITE  (6,7220)  (NBN(l) , 1*1 , MB) 

WRITE  (6,7225)  N8YMB, N8YM6, NPLff , NBNO, NDERV 
WRITE  (6,7150)  NNOUT 

IF(NNOUT.CT.O)  WRITE (6, 7 180)  (MOUT ( I ) , MOOT ( I ) , 1*1 , NNOUT) 
WRITE  (6,7105)  NKP 

IF  (NKP.CT.O)  WRITE  (6,7160)  (KPG ( I) , KPB ( I ) , I»1 , NKP) 

WRITE  (6,11600)  NL,XLP 
IF(NPtT.EQ.O)  WRITE(6, 7230)  THETAO 
IF(NPLIT,EQ,1)  WRITE  (6,  7260)  THETAO, A 
IF  (N0ERV.EC.2)  CO  TO  1180 
00  1160  I«1,NL 

1160  WRITE  (6,11700)  I , HM ( I ) , RHOM ( I ) , EX ( I ) , ET ( I) , XXNU(i) , THNU ( I) , 
1 CXT(I),8AT(I),8AC(I) 

IF  (KTVPE,NE.J,AND.KTYPE.NE,4)  CO  TO  1190 
IF  (KOAM.NE.O)  CO  TO  1190 
WRITE  (6,11900)  EC,OC,DC 
CO  TO  1190 

1180  WRITE  (6,7260)  (HM(I ) , RHOM ( I ) ,EM(I),I*1,NL) 

WRITE  (6,7300)  TNU,SIG0,EP,EP8IF 
1190  WRITE  (6,7400)  ( (FC( I , J) , 1*1 , MB) , J*1,MC) 
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WRITE  (6,  WO)  OELTIM*  T3T  OP 'PRINT 
C 

2100  I»0 

MBMBaO 

DO  2110  Mal#M6 
NX  • MBM(M) 

DO  2150  Nal'MB 
NN  • NBN(N) 

MUBC(N#M)al 

IPU.tO.NNOUT)  BO  TO  2150 
DO  2110  Jat'NNOUT 

IF(MM,EB.MOUT(J).ANO,NN.EO.NOUT(J))  SO  TO  2120 
2110  CONTINUE 

BO  TO  21  SO 
2120  MUSE(N,M)«0 
!•!♦! 

•0  TO  2110 
21  JO  M6MBNM8MB*! 

2150  CONTINUE 

MBNIiWMBMB 
00  2100  Mat 'MB 
MM  a MBM(M) 

CC1(M)  a MM 
2200  CCS(Mi)  a MM  ♦ 1 
KJaiBO.O/TMiTAO 
IF  (NPLIT.EB.O)  XJaPI/THETAO 
00  2300  Nal'MB 
NN  • NBN(N) 

CC2 (N)  a NN 
2300  CCO(N)  a NN  ♦ 1 
RETURN 


C 

7000  FORMAT (BI12) 

7100  FORMAT (BF12« 1 ) 

7150  FORMAT  (9HONNOUT  a 15) 

7170  FORMAT  (2AH1 INPUT  DATA  FOR  DEPROP  ) 

7100  FORMAT  (2X«) 

7105  FORMAT  (TMONKP  a 13) 

7200  FORMAT ( 10H0M8 

12/10H  MBAR  a 12/1  OH  NBAR  a I2/10H  LB  AR 


a 12/ 1 OH  MB 
a 12) 
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r<tut  FUWMAI  (10H0HGM  ■ ( 
7220  FORMAT  C l OHONBN  ■ ( 
/ 225  FORMAT  (1QM0N8YMG  ■ I 
3 10H  M8ND  ■ 12/ 

7230  FORMAT (17H  THE T AO#  IN 
7260  FORMAT ( 1 7H  THETAQ,  DE6 


(1015)) 

(1015)) 

I2/10H  M8YMB  ■ 12/iOHONPLT 
10H  NOERV  • 12) 

» E16.8) 

C * E16.S/17H  A # IN 


7280  FORMAT  (12N0  HM#  I N , AX , 21 HRMQM,  LB«8EC**2/ 
1 7 HEM#  PSI/ ( JE1 7 .8) ) 

7300  FORMAT  U7H0TNU  « E16.B/17H  8100#  P5I 

1 P81  * E18.6/17H  EP8IF#  IN/IN  « E16.6) 

7400  FORMAT(6H0FC  « / (5214,6) ) 

8200  FORMAT (15HOOELTIM,  SEC  • E16.6/1SH  TSTOP,  SEC 
1 ■ £16,8) 

1U00  FORMAT  UOHONl  « I2/17H0XLP,  IN  • El 


• E16.8) 


LB«9EC**2/IN**4#4X, 

8100#  PSi  « E46.6/17H  EP, 
« £16,8) 

STOP#  SEC  ■ E16.S/15H  PRINT 


11700  1 

FORMAT 

(6H0LAYERI3/27H  HM# 

IN 

• 

1 

27H 

RHOM#  LB-SEC**2/IN**4 

a 

£16,6/ 

2 

27H 

EX#  PSI 

a 

£16.6/ 

3 

27H 

ET,  PSI 

a 

E16.S/ 

4 

27  H 

XXNU 

a 

£16,6/ 

5 

27M 

THRU 

a 

£16.6/ 

6 

27H 

6XT#  PSI 

a 

£16,6/ 

7 

27  H 

SAT#  PSI 

a 

£16.8/ 

8 

27M 

SAC#  PSI 

a 

£16. 6) 

11900  1 

format 

(11H0EC#  PSI  • £16,6/1 

1H 

GC,  PSI  « £16. 6/ 

1 

1 1 H 

OC#  IN  • E16.6) 

El  6,8) 

* £16,8/ 
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•DECK  D8ET2 

SUBROUTINE  0SET2 
C 

•CALL  CBLK1 
•CALL  CSLK2 
•CALL  CBLK3 
•CALL  CBLA« 

•CALL  CBLK9 
•CALL  CBLK7 
•CALL  CBLKB 
•CALL  CBLKA 
•CALL  CBLK10 
•CALL  CBLK1I 
•CALL  CBLK1S 
•CALL  CNOVA 
•CALL  CBLANK 
C 

IF  (NOKRV.EQ.l)  60  TO  2710 

call  lcsind 

IF  (LBAft'f Qt 0)  60  TO  3000 
IF  (RTYFE.NE.S)  60  TO  2500 

C SET  UF  EQUIVALENT  LAYER  FOR  HONEYCOMB  METAL#  MDERV  * 2. 
NELF  » 2 
H2«HM(J)*0,5 
H01«HM(3)*HM(2) 

HDitHM(l) 

ZB(l)«*Ht«0.5*H02 
ZB(t)  » H2  • ,5*H01 
SUNRHSO.O 
HORO.O 

00  2400  1*1 » 3 
H1*HM(1) 

SUMRH«8UMRHtRHOM(I)*(Hl-HO) 

2400  H04HI 
NL*I 

EQCRHM (3) •HM(2)»HM(l) 

EQH»(HM(5)»MM12)»HM(1) )«80RT (3,0*HM(l)*(HMC3)«HMC2)))/EQC 

RH0»8UMRM/EQH 

H*|QH 

HM( 1 ) *M 
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non 


1 11  11  — wHiouan  1 1 ui  ■■■»»  ■ .. . ...t  ji.jWMmiiuv,-.tiWpiuii..w  1 JIW- 

1 


EL*EM(l}*EQC/ESlH 
£P*EP*EBC/EOH 
S160a9160*EQC/EUH 
EM( 1 ) aEL 

IF(AQAM,E6.2)  60  TO  1705 

TC«UT(imPSIP 

60  TO  2705 

C SINGLE  LAYER,  NDEWV  * 2, 

2500  H * HM C 1 ) 

EL»EM  C 1 ) 

RHOaRHOM(l) 

H2«0.5*H 
ZB ( 1 )t*H2 
ZB ( 2)  ■ H2 
NELP  * 2 

IF(AOAM-l)  2600,2700,2705 
C NO  0AMA6E 
2b00  TCRIT(1)bSIG0 
CCRIT (l)a8I60 
NELR  • 1 
GO  TO  2705 
2700  TCRIT(i) aEPSIF 

IFUTYPE.EO.l)  60  TO  2705 
TCRIT(I)  * SI60 
NELP  • I 

2705  OT  a OELTIM 

F3«RMO*U,0  • TNU**2)/EL 
Ft  • Ht*l/(i2,0*F3) 

F 2a  8I60/RH0 

CALL-  OTSTEP  (Fl,F2, F3,F3,F1) 
IF  (DT.GT.0.0)  OELTIM  a OT 
60  TO  2750 

NOERV  a I, 

COMPOTE  HBAR, 

2710  A 1 a 0. 

A2  ■ 0, 

AS  • 0. 

A«  a 0. 

AS  • 0. 
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1 


! 


AO  • 0, 

A7  ft  0, 

Aft  ft  0. 

HO  ft  0. 

00  if  IS  l«l,NL 
HI  ft  HH(U 

• 2i  » l./(l,  • XXNU(I)*THNUU)) 

01 t • EX(I)*»22 

ftii  • et (D*fti2 
IK  • XXNU(  I ) *822 
•IS  • OXT(S) 

• XL  ( I ) ft  BU 

STUD  • B2{ 

01  ft  HI  • HO 

02  • Ml**2  • H0**2 

A1  ft  A1  ♦ 111*02 

A2  • A2  ♦ BllftOl 

A3  • AS  ♦ 022*02 

Aft  ft  Aft  * 022*01 

AS  ■ A5  ♦ 013*02 

At  a AO  ♦ 033*01 

Af  ft  Af  ♦ 012*02 

AO  ft  AO  * 012*01 

2715  HO  • HI 

HO l l ft  •S*A1/A2 
H022  ft  *5*03/04 
H8I3  • *5*A5/A6 
HO  1 2 ft  *S*A7/A0 

HOAR  ft  *23*(H0n  ♦ HB22  ♦ HB33  ♦ H012) 

•HOBO  • 0, 

CHI  1 ft  0, 

CHI 2 ft  0, 

CM22  « 0* 

CM33  • 0. 

• MU  ■ 0* 

• Ml 2 ■ 0* 

•M22  ft  0, 

•M33  • 0* 

OMll  ft  0. 
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DM 1 ? * o, 

DM22  a 0. 

DM33  ■ 0. 

FaO.O 
HO  ■ 0. 

DO  2720  I«1,NL 
HI  s HM II) 

BTT  « BTL(l) 

KNOB  « XXNU( I ) *BT  T 
BXX  » BXL(I) 

CXTL  a 6XT(I) 

«HOL  ■ WHOM ( I ) 

H 1 1 a HI  • HO 
F*  F ♦ SAT ( I ) *H1 1 
CM J 1 c CMU  ♦ BXXftHl  1 

CM12  a CM12  ♦ XNUB*H1 1 

CM22  a CM22  t BTTaHU 

CM33  a CM33  ♦ 6XTL*HU 

RHOBR  a RHOBR  ♦ RH0L*H11 
HI  2 a Hl**2  • H0**2 
HI  2D  a H12  • 2,*HBAR*HU 
FMli  a FM11  ♦ BXX*H12D 

FM1  2 a FM12  ♦ XNUBM12D 

FM22  a FM22  ♦ BTT *H12D 

FM33  a FM33  ♦ 8XTL*H12D 

HI 3 a Hl**3  • H0**3 
HI 30  a HI 3 - 3,*HBAR*H12  ♦ 
DM11  8 DM11  ♦ BXX«H1 3D 

DM12  a DM12  ♦ XNU8*H 1 30 

DM22  a DM22  ♦ BTToHIJO 

DM33  a DM33  ♦ 6XTL*H130 

2720  HO  a HI 

0H3  a l./HMCNU 
OA  a l./A 

02A2  • 1 •/ (2a*A**2) 

03A3  a 1,/(3,«A*«3J 
CMU  • CMU*OA 
CM12  a CMltaOA 
CM22  a CM22*0A 
CM33  a CM33*0A 


3V*HBAR**2*H11 
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FM11  « FM1U02A2 
MU  ■ PM12*02A2 
FM 22  « FM22*02A2 
FHSS  • FM31*02A2 
0M11  ■ DMU*03A3 
ONI  * ■ 0M12-03A3 
0*22  ■ 0M22*03A3 
0M33  ■ 0M3J*O3A3 
RHQBR  ■ RHQBR*0H3 
RHO  • RHOBR 
H • MN(NL) 

F3»RH0/0H3 

F1«0M22*A**3/F3 

F2*F/F3 

F4  « F3/(A*CHU) 

F5  a 0*1 1*A**3/F3 
F3»FS/(A*CM22) 

OT  a DELUM 

CALL  QT8TEP  (F1,F2,F3,F4,F5) 

IF  (9T.6T , 0,0)  DELTIM  a OT 
NELP  • 1 
NZP  * 2 

IF  (KTtPE.EO.5)  60  TO  2730 

NLZU)  a 1 

ZCU)  » -HBAR 

NLZ (2)  » 1 

ZC (2)  • ZCC1)  ♦ M 

IF  UTYPE.LT, 3)  60  TO  2745 

NLZ (2)  a 3 

ZCU)  • ZCU)  ♦ , 5*HM U ) 

ZCU)  > ZC  (2)  • ,5*(HM(3)  • HM(2) ) 
60  TO  2749 
2730  HT  a •NBA* 

HS  * M*U) 

00  2740  latfNL 
NLZtiftl  • 1)  * I 
NLZ (2*1 ) • I 
ZC (2*1  • n » HT 

IF  II.6T.1)  M8  * HM(I)  • HM(I-l) 

HT  i HT  « H9 


CO  TO  2600 


2/ 40 

ZC(2*D 

> HT 

NZP  a 2 * 

NL 

2/45 

IF  UOAM 

.EU, 2i 

00  2747 

i»i»nl 

TCHITd) 

■ SAT 

2747 

ccrituj 

» SAC 

2750 

CONTINUE 

2600 

RETURN 

3000 

KERR  ■ I 

RETURN 

END 
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• DECK 

D8CTS 

p 

SUBROUTINE  DSET5 

•call 

CSLK1 

•CALL 

CSLXC 

•call 

C6LKS 

• CALL 

CSLK4 

•call 

C6LK5 

•CALL 

C6LK7 

• CALL 

CSLKO 

•call 

C6LK9 

• CALL 

C6LK10 

•CALL 

C6LK11 

• CALL 

CBLK1J 

• CALL 

CBLK14 

• CALL 

CNOVA 

• CALL 

C6LANK 

c printout  description  of  deprop  data 
c 

2800  MRITC<6,9J00) 

IF(NPlfT.EQ'O)  NRITC«6,9400) 

IF(NPw/T.E8,U  MRI T£<6,  9500) 

SO  TO  (2020,2640,2660,2660,2900) , KTTPE 
2620  MRITC  (6,9600) 

60  TO  2950 
2640  MRITC  (6,9700) 

60  TO  2950 
2660  MRITC  (6,9600) 

60  TO  2950 
2860  MRITC  (6,9620) 

60  TO  2950 
2900  MRITC  (6,9640) 

2950  IF  (N6NO.E6.1cOR,NSN0.CO,3,OR,NBND,EO,6) 
IF  (N6NDVC0,2,0R.NSND,CS.4,0R,NBN0.C0.6) 
IF  (N6N0,Ce.5,0R,NBND,EQ.7,0R.NBND.CQ,9) 
IF  ( NBNO. E9. 1, OR. NSND.EQ. 4. OR. NBNO.EQ. 5) 
IF  (M6N0,C8,2,OR,NSNO.EO,3,OR,N8NO,EU.n 
IF  IN6NO,EQ.6,QR,NSNO.CQ.6.OR,N0ND,CQ.9) 
IF(NOCRV.eO.l)  MRITC(6, 10(00) 


MRITC  (6,9900) 
MRITE  (6,9920) 
MRITC  (6,9940) 
MRITE  (6,9960) 
MRITC  (6,9960) 
MRITE  (6,10000) 
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IF (*DfcWv.EQ.2)  WRITE (b, 10200) 

WRITE (6, 10000)  M6,MB,MBAR,NBAR,LBAR 
WRITE  (6' 1 0620) 

00  2970  M*1,MG 
MM  a MGM(M) 

00  2970  N*l,MB 
NN  a NBN(N) 

IF  (MU8E(N,M),EQ,C)  GO  TO  2970 
WRITE  (0,10630)  MM, NN 

2970  CONTINUE 

WRITE  (6,10650)  XL? 

1F(NPLT«EQ.0)  WRITE(6, 10900)  THETAO 
IF(NPL‘T.E0,1)  WRITE (0, 11000)  THETAO,  A 
IF  (N0ERV.EQ.2)  GO  TO  3500 
WRITE  (6,12050)  HBAR,  ( I , Ial,  NLi) 

WRITE  (6,12100)  (HM(I),Ial,NU 
WRITE  (6,12200)  (RHOM ( I ) , Ial , NL) 

WRITE  (6,12300)  (EX(I) , I*l,NLt) 

WRITE  (6,12*00)  (ET(l),l8l, NLO 
WRITE  (6,12500)  (XXNU(l) , Ial ,NL) 

WRITE  (6,12600)  (THNU( I) , 1st , NL) 

WRITE  (6,12650)  (6XT(I),lal,NL0 
IF  (KTVPE.NE.1.AND.KTVPE.NE.3)  GO  TO  3300 
WRITE  (6,12900)  (BAT (I)  , 1*1,  NLt) 

WHITE  (6,13000)  (SAC (I ) , 1*1 , NL) 

GO  TO  3*00 

3300  WRITE  (6,12700)  (BAT (I ) , 1*1 , NLO 
WRITE  (6,12600)  (SAC ( I ) , 1*1, NLO 

3*00  IF  (KTTPE,NE,3.AND,KTYP£,NE,4)  GO  TO  3600 
IF  (KOAM.NE.Q)  GO  TO  3600 
WRITE  (6,13100)  EC,GC,OC 
GO  TO  3600 

3500  WRITE  (6,11100)  H, RHQ, EL, TNU, SIGO, EP, EPSIF 

3600  WRITE  (6,11200)  ( (FG (N, M) , N*1 , MB) , M*1 , MG) 
WRITE (6, 11300)  OELTIM'TSTOP, PRINT 
IF  (NOERV.20,1)  60  TO  4020 
DO  4000  K*1 ,LBAR 
ZH(K)*BX(K)*H*0,5 
If (K)*ZH(K)/A 
Z6(K)*6X(K)a*2 
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4000  CONTINUE 

ZA(1)  • ZB(1)/A 
2 A 12)  • ZB(2)/A 
1NZU)  • 1 
INZ (!)  • LBAR 
4020  NNBYMO  « 0 
NNBTMB  • 0 

IF  (N8ND,E8.5,0R.NBND.£Q.7.QR.NBND.Ea,9)  NNBTNB  * 1 
IF  (NBND.E0.6.0R.NBND.6E.6)  NNBTMB  ■ 1 
IF  (NNBYMB.tO.l.AND.NSYMB.EO.O)  WRITE  16# 1 3200 ) 

IF  (NN8YMB.E0.1.AND.N8YMB.E0.Q)  WHITE  16,13200) 

NBT  * MBAR 
NBT  ■ NBAR 
NB  • NBT 
NB  • NBT 

IF  (NBVMB.CO.l)  NB  « (N8Wl)/2 
IF  (NBYMB.EO.l)  NB  a (NBTtl)/2 
NY2  a J*M6MB 

F1M  a PI/FLOAT (2* ( MB AR«1)) 

PIN  a PI /FLOAT (2*(NBAR*1)) 

IF  (NBYMB.EO.l)  PIM  ■ 2,*PIM 
IF  (NBYMB.EO.l)  PIN  • 2.«PIN 
RaA/H 

XLaXUP/(PI*A) 

XLial.O/XL 

XL2aXUM2 

XLSat.OaXLl 

XL«at«B*XL2 

XLSa*L»4B 

XL?*1«0/XL2 

CNB  a XL**2 

CN9  a CNB/ (2.0*R) 

IF  (NOCRV.EQ.l)  BO  TO  4040 
Cl  a 1.0/ (Aa*|*XL**S) 

C2  a (XJ/A)**2/XL 
C3  a (X J/A) **2*X J 
C4  a XJ/(A*XL)**2 
C9  a •H**2/2.0 
Cfc  a H**>/0.0 
OCLR  • PIM*XLP/PI 
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4040  CN10  • 2,*CN6*R 
CNU  • CN10 
Cl  a D*U*A/XL**3 
C2  * C0M12  ♦ 4.0*DM33)*A*XJ**2/XL 
C 3 « DM22*A*X J**3 
C4  a (DM12  ♦ 4.0*DM33) *A*XJ/XL**2 
C5  a •4,0*OM33*XJ*A**2/XL 
4QbO  XJ2  a X J**2 
XJ3aXJ«XLl 
XJ4a2,o*XJ 
XJ5a2.0*XJ3 
OPRTlaPRINT*OELTIM 
XLPlaXLS 
XLP2«2.0*XL*P 
PPL1  a 1,/(PH0*XLP**2) 

C SIMPSON  S RULE. 

C MBAR  AND  NBAR  MUST  BE  ODD  NUMBERS  POP  POLL  PANEL* 
00  4100  lal.MBAR 
P a I-l 

CAM(I)  a FaPlM 

XG (I ) a CAM(I)*XLP/P1 

PIMA(I)  a PRL*4.«PlM/3. 

IP  (N8YMG.E0.1)  PIMA(I)  a ,5*PIMAU) 

IF  (UM)/2,EU.I/2)  PIMA  (1)  a 2,*PIMA(I) 

4100  CONTINUE 

PIMA(l)  a PIMA  Cl)*. 5 
PIMA (MBAR)  a PIMA (M8AH) * ,5 
DO  4200  lal.NBAR 
P a 1*1 

BETR(I)  a PaPIN 

XB ( I ) a BETR(I)*THETAO/PI 

PINAU)  a 4,*PIN/3, 

IP  (NfVMB.EQ.l)  PINAU)  a ,S*PINA(I) 

IP  ( (1*1)  /2.E0.I/2)  PINAU)  a PINA ( I ) *2. 

4200  CONTINUE 
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pinacd  • .9o*piNA(n 

PINA (MB AM)  • PINA (NBAR) *«5 
00  4400  Ial,MBAR 
00  4490  Jal,NBAR 
NBUBtt Jt I)  • 0 
4400  NUBEIJ'l)  • 2 
II  a 0 

00  4430  Ial,N6T 
00  4410  J*1 # NBT 
IP  II2«E9*NKP)  80  TO  44)0 
00  4410  Kal,NKP 

IP  (I.EQ.KPS(X),ANO.J.EO.KPB(K) ) 60  TO  4420 
4410  CONTINUE 

00  TO  44J0 
4420  NUSCUtl)  • 3 
II  • II  ♦ 1 
4430  CONTINUE 

Nt  • NBAR  • NtTMB 
K a 0 

00  4440  jBtfN! 

K •.  K ♦ I 

IP  (NUtE( J«  t) •0T«2)  NBU8E ( J,  1 ) a *K 

IP  (N«YM6,E9tl.AN0,NUUWf«BAR).0T,2)  NBU3E ( J#  MBAR)  * .K-N1M 
4440  CONTINUE 

N1  a MBAR  • NOT MS 
Kao 

00  4 450  Ia2f N1 
K a K ♦ 1 

IP  (NUBE(1,I),6T.2)  NBU8E (1*1)  a K 

IP  (NBYMB,EQ.1.AN0.NUSE(NBAR,I),6T,2)  NBUSEtNBAR,  I)  a KtNfl 
4450  CONTINUE 

NBUSUi'l)  a 101 

IP  (NBYMB.EO.l)  NBUBE(NBAR,1)  a 102 
IP  (NBVM6,C0,1)  NBUSE(l«MEftR)  a 103 
IP  (NSVMBaNiVNB.EQ.I)  NBUBE (NBAR, MBAR)  a 104 
NRC  a l ♦ NSYMB  ♦ 2»NBYM9 
Ifao 

00  4500  Mai, MS 
XI  a MBM(M)  ♦ 1 
00  4500  Iai,N6T 
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u*im 

»1H6UI)*81NU1*6AMU) ) 
C0Se(maC09(Xl*6AMU)) 
»lN26in)«SlN((Xl«l,0)*6AM(I)  j 
C0926<(II)*C03((X1»i,o)*&AM(I)) 

4500  CONTINUE 
II«0 

00  AbOO  N»1,MB 
XI  c NBN(N)  ♦ 1 
00  AbOO  J«1 > NBT 
XI»II»1 

SINB(U)*9IN(Xl*BETft(J)) 

C08B(XI)*C03(X1*BETR(JJ) 

81N2B(U)b8IN((X1»1,0)*BETR(J)) 

C082B  ( I I)«C03((Xl»|,0)*BETR(J) ) 

AbOO  CONTINUE 
C 

CALL  BOLT 
C 

ABOO  K»0 

OAAl.O/A 
00  5200  1*1 f N&T 
00  5200  Jal'NBT 
IF<NU8EU,I).EU.0>  90  TO  5200 
KAXM 

OmO(K)aO«0 
0*9U)«Q,0 
DAB (K J *0«  0 
00  5100  M«lrH9 

1 ) *N6T  ♦ I 
00  5100  N«1,MB 

IF (MUSE (N,M)(EU,0)  90  TO  5100 
NN«(NM)*NBT  ♦ J 
FG«N«F9(N,M)*0A 

OMO(K)  • F6MN*FFl(MN)*FF5(NN)  ♦ DNO(K) 
OHft(K)  • F9MN*FP2 (MN) »FF5 CNN)  ♦ DNB(K) 
D*B(KJ  A FSMN*FPl(NM)*Ffb(NN)  ♦ DNBOO 
5100  CONTINUE 
52o0  CONTINUE 
N6NBT  ■ K 
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NONA  ■ N6»NB 
RETURN 

9300  FORMAT (1*1,29*, IJHD  ERROR  /15H0PANEL  ANALYZEO) 
9400  FORMAT (7M  FUAT) 

9900  FORMATION  CURVED) 

9600  format (tin  metal*  single  LATER) 

9700  FORMAT (26H  PLASTIC*  SINSLC  LATER) 

9600  FORMAT ( 1 9H  METAL*  HONETCOMS) 

9S20  FORMAT  (I1H  PLASTIC*  HONETCOMS) 

9640  FORMAT  (22H  PLASTIC*  MULTILATER) 

9900  FORMAT  (37H  CLAMPCO  • CLAMPEO*  GAMMA  DIRECTION) 

9920  FORMAT  (SSH  SIMPLE  • SIMPLE*  GAMMA  DIRECTION) 

9940  FORMAT  (SEN  CLAMPCO  • SIMPLE*  GAMMA  DIRECTION) 

9960  FORMAT  (16H  CLAMPEO  • CLAMPEO*  BETA  DIRECTION) 

9960  FORMAT  CS4M  SIMPLE  • SIMPLE*  SETA  DIRECTION) 

10000  FORMAT  (SSH  CLAMPEO  • SIMPLE*  BETA  DIRECTION) 

10100  FORMAT (26H0RESP0NSC  OPTION  . ELASTIC) 

10200  FORMAT (34H0RESP0NSE  OPTION  • ELASTIC-PLASTIC) 

10600  FORMAT ( 1 7H0S TRUC TuRAL  MOOEL/ 

1 4 7 H NUMBER  OF  GAMMA  MODES  (MG)  ■ 13/ 

2 47H  NUMSER  OF  BETA  MOOES  (MB)  ■ 13/ 

3 «7M0  NUMBER  OF  GAMMA  INTEGRATION  P0INT9  (MBAR)  • 13/ 

4 47H  NUMBER  OF  SETA  INTEGRATION  POINTS  (NBAR)  • 13/ 

5 47H  NUMBER  OF  Z INTEGRATION  POINTS  (LBAR)  * 13) 

10620  FORMAT  (24HOMODAL  COMBINATIONS  USED) 

10B30  FORMAT  (3**214) 

10650  FORMAT  (19N0  LENGTH  OF  PANEL*  IN  (XLP)  » E16.6) 

10900  FORMAT ( J9H  MlDTH  OF  PANEL,  IN  (THETAO)  « £16. 6) 

11000  FORMAT (ISH  8UBTEN0E0  ANGLE*  OCG  (THETAO)  ■ £16.8/ 

2  SSH  RAOIUS*  IN  (A)  s E16.6) 

11100  FORMAT ( J5N  THICKNESS*  IN  4 £16.8/ 

1 SSH  DENSITY*  LS-8EC*«2/IN*«4  * £16.6/ 

2 SSH  ELASTIC  M00ULU8*  PSI  * E16.6/ 

3 SSH  POISSON'S  RATIO  ■ E16.S/ 

4 SSH  TIELO  STRESS*  PSI  • E16.S/ 

5 SSH  STRAIN  HANOCNING  SLOPE*  PSI  « 116.8/ 

6 SSH  ULTIMATE  STRAIN*  IN/IN  (EPSIP)*  E16.8) 

11200  FORMAT (26H0INITIAL  IMPERFECTIONS*  1N/(5E14.6)) 

11300  PORMAT(17HOTIM£  INFORMATION/ 
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1 «2H 

2 42M 

3 42H 

12050  FORMA  I 


INTE6RATIQN  STEF  SIZE,  SEC  (OELTIM) 
STOP  TIME,  SEC  (TSTOP) 

PRINT  FREUUENCY  (PRINT) 
(40HOCOOROINATE  SURFACE  POSITION 


1 13H0L.AVER 
12)00  FORMAT  127H 
12200  FORMAT  (32H 
12300  FORMAT  (33H 
12400  FORMAT  (40H 
12500  FORMAT  (22H 
12600  FORMAT  (26H 
12700  FORMAT  (JIM 
1 2600  FORMAT  (3SH 
12900  FORMAT  (26H 
15000  FORMAT  (52M 
12650  FORMAT  (21H 
13100 


« £16.6/ 

« 116, a/ 

• E16,S) 
(HSAR),  IN 


£16.6/ 


FORMAT 

FORMAT  (62HOCORE 
1,  PSI  « €15.6/ 

2 63H  SHEAR  MODULUS 

3 £14.6/ 

4 62H  CORE  CELL  SIZE 

5 E15.6) 

13200  FORMAT  (41H0**  MARNIN6 
ENO 


NUMttER, 22X , 41 1 5/ ( 3) X , 41 1 5) ) 

CUMULATIVE  THICKNESS,  IN, 13X,6El5.6) 

MASS  DENSITY,  LB*SEC««2/IN**4,6X,6E15.6) 

MODULUS  OF  ELASTICITY  • X,  PSI, 7X,6E15,6) 
MODULUS  OF  ELASTICITY  • THETA,  PSI  6E15.6) 
POISSON'S  RATIO  • X, 18X, 6E15.6) 

POISSON'S  RATIO  - THETA, 14X,6ElS,6) 

TENSILE  ULTIMATE  STRESS,  PSI, 9X,6E15,6) 
COMPRESSIVE  ULTIMATE  STRESS,  PSI,5X,6E15,b) 
TENSILE  YIELD  STRESS,  PSI, 12X, 6E15.6) 
COMPRESSIVE  V1IL0  STRESS,  PSI , SX, 6E1S.6) 

SHEAR  MODULUS,  PSI, 19X,6E13,6) 

MODULUS  OF  ELASTICITY  PARALLEL  TO  CORE  DEPTH 


(EC) 


OF  CORE  (60,  PSI 
(DC) , IN 


•*  INCONSISTENCY  IN  SYMMETRY) 
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•DEC*  DT9TEP 

SUBROUTINE  DTSTEP  (FI  * F2, F Jr F4, F5) 

INIS  SUBROUTINE  COMPUTES  AN  APPROXIMATE*  CONSERVATIVE 
TIME  STEP  FOR  DEPROP. 

NOTE!  IT  IS  ASSUMEO  THAT  FOR  A FLAT  PANEL  T HE T AO  ,LT.  XL. 


•CALL  CBLK1 
•CALL  CBLK2 
•CALL  CBLK13 
•CALL  CNOVA 
C 


c 


BN  ■ NBAR 

IF  (NSTMB.EQ.O)  BN  • 2*NBAR  - I 
THITR  ■ PHTHETAO/ISO, 

BM  ■ MBAR 


CMaO.O 

IF (NBNO  *EO. I 

IF(NBN0.EQ,5 

CN«0.0 

IF (NBNO.EQ. 1 


DTI  X 
0T2X 
0T3X 
OT  AX 
0T5X 


0)  BM  a 2AMBAR  • 1 

.OR. 

MONO. IQ. 3 .OR. 

NBNO.EQ. S) 

CMaO . 30 

.OR. 

NBNO.EQ. 7 .OH. 

NBNO.EQ. 9) 

CMaO.15 

.OR. 

NBN0.EQ.4  .OR. 

NBNO.EQ. 9) 

CNaQ , 30 

.OR. 

NBNO.CE.S)  CNaO 

.15 

l.OES 

1.0E6 

l.OES 

l.OES 

l.OES 


IF  (NPLT.EO.O)  DT2X  a 
IF  (NPLT.EQ.l)  0T3X  • 
IF  (NPLT.EQ.l)  DT AX  a 


check  all  modal  combinations, 
DO  200  Mai, MG 
BARM  * MSM(M) 

00  200  Nai , MB 

IF  (MUSE(N.M) ,EQ.O)  GO  TO  200 
BARN  * NBN(N) 

0 ARMXaB ARM  ♦ CM 


THE T AO* SORT (P3)/(BN  • 1.0) 
A*THETR*8QRT(F3)/(0N  • 1.0) 
XLP*80RT(F«)/(BM  • 1,0) 


BARNX  a BARN  ♦ CN 
IF(NPLT.EU.l)  so  to  too 


DTSTEP 


c 


FLAT  PANEL 

IL*N  ■ (BARMX*PI/XLP)**2  ♦ (BARNXaPI/THETA0)**2 
Oil  a PI  /S8RT(ELMN*(F2  ♦ ELMN*F 1 ) ) 

Oil  • DTI  / 25.0 
IF  (DTI. LT. DTIX)  DTIX  ■ 0T1 
GO  TO  200 
C 

C CURVED  PANEL 

100  EL«IN  ■ BARMX«P1*A/XLP 
EKMN  ■ BARNX*P1/THETR 
DUM  • (ELMN**2  ♦ EKMN**2)**2 

DTI  a P I • A «S8RT(P3/(Q.5*(1.0  ♦ ELMN**2)  • 0.5*S8HTU1.0  • 

1 ELMN*«2)**2  t «.0*(0.30«ELMN)**2))J 
DTI  a DU/35.0 

DTE  a P I a A /SORT (FI*DUH/A**2  ♦ ELMN«««/ (F3*DUM) ) 

DTE  a 0T2/35.Q 

ELMN  a S8RT(0UM)/A**2 

DT5  a PI  /$QRT(ELMN«(F2  ♦ ELMN*F5) ) 

0T5  a DT5/2S.0 
IF  (DTI.LT.0T1X)  DUX  a 0T1 

IF  (DT2.LT.DT2X)  0T2X  a 0T2 

IF  (DTS.LT.0T5X)  0T5X  a DT5 

200  CONTINUE 

IF  (NPLT.EQ.O)  DEL  UM  a AMIN1  (DT  1 X,  DT2X) 

IF  (NPLT.EO.l)  DELTIM  a AHIN1 (DT 1 X , DT2X, DTIX , DT 4X , 0T5X) 

IF  (NOBUG.  GT.O)  WRITE  (b,1000)  DT  1 X , DT2X , 0T3X , DUX  , 0T5X 
RETURN 
C 

1000  FORMAT  (30H0DEPRUP  TIME  STEP  CALCULAT I0NS/5E1S.6) 

END 
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•DECK  HIM 

SU9R0UT  INE  HIM  (*,  NfQ,  OELT  IM,  TIME,  VXO,  X4,  AX«) 

SPECIAL  INTEGRATION  METHOD  FOR  2ND  ORDER  DIFFERENTIAL'  EUUATIUN3 
PHICH  HAVE  NO  DAMPING.  CENTRAL  DIFFERENCE  SCHEME. 

COMMON/CHIM/  X3(147),DTSQ 
DIMENSION  X4(l),AX4(l),VXQll) 

c 

IF  IX.GT,  1)  GO  TO  200 
DTSB  a DELTIM**2 
K a 2 

DO  100  iai'NEO 

100  XJU)  • X4(I)  • VXO(l)«DELTIH  ♦ 0 . 5*D  TSG*  AX4 ( I ) 

C 

200  DO  300  1st  , NEW 

X • 2, 4X4(1)  • X 3 ( I ) ♦ OT80* AX4  (I) 

XI(I)  a X4CI) 

300  X4(I)  8 X 

TIME  a TIME  ♦ DELTIM 

RETURN 

END 


ouou  t 


• UfcC*  LEStNO 


SUBROUTINE  LE6EN0 

THIS  ROUTINE  CONTAINS  THE  LE6ENDRE  ZEROES  A NO 
NEIfiHUNS  FACTORS  FOR  18AR  .Li.  IS. 

CALL  CBLK1 
•CALL  CBLK3 

DIMENSION  CCX(7,14),  CH0(7,14) 

C 

DATA  CCX/ 7 *0.0, 0.57 7 350269 16963, 6* 0.0, 0,7 7459666924 148, 

1  6*0.0, 0.86113631159405,0. 339961043*6466, 9*0.0, 0.90 617 95 A *9 1666, 
i 0.53646931010966,9*0,0,0.93246991420315,0.66120935646626, 

3 0.23861916608320,4*0.0,0.94910791234276,0.74153116999939, 

4 0.40564515137740,4*0.0,0.96026965649734,0,79666647741363, 

5 0.52553240991633,0.16343464249565,3*0.0,0,96616023950763, 

6 0.636031 10732664, 0.61 337 143270059, 0.32429342340361, 3*0.0, 

7 0,97390652651717,0.66506336660696,0.67940996629902, 

8 0.433395394 12925,0. 14887433696163, 2*0, 0,0. 97622665614606, 

9 0.66706259976609,0.73019200597405,0.51909612920661, 

1 0.26954315595234,2*0.0,0.96156063424672,0.90411725637047, 

2 0.76990267419430,0.56731795426662,0.36763149699616, 

3 0.12523340651147,0.0,0.96416305471859,0.91759639922296, 

4 0.80157809073331,0.64234933944034,0.44649275103645, 

5 0.23045631595514,0.0, 0.96626380669661 , 0*92643466366357 , 

6 0.82720131506976,0.66729290461166,0 . SI 524663635815, 

7 0.31911236692789,0.10605494670734/ 

C 

DATA  CHO/ 7 *0.0, 1.0,6*0.0,0.55555555555556,0.88866666866689, 

1 5*0.0,0.34765484513745,0.65214515486255,5*0.0,0.23692686505619, 

2 0.47662667049937,0.56866666868669,4*0.0,0.17132449237917, 

3 0.36076157304614,0.46791393457269,4*0.0,0.12946496616887, 

4 0.27970539146926,0.36163005050512,0.41795916367347,3*0.0, 

5 0.10122853629038,0.22230103445337,0.31370664567769, 

6 0.36268378337636, 3*0. 0,0. 61 27438836 1574C»1, 0.1606461 6069486, 

7 0.26061069640291,0.31234707704000,0.33023935500126,2*0.0, 

8 0.666 7 1344*06666E«1, 0.149451 34915058, 0,21 906636251 596, 

9 0, 2692667 1931 000, 0.2955242247 1 475, 2*0 .0,0. 55666567 1 161 74E«1, 

1 0.12558036946491,0.16629021092773,0.23319376459199, 

2 0.26260454451025, 0.27  29250667 7790,0.0,0.47 1 7 531636651 2E*1 , 
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3 0. 1 0693932599532, 0.16007 e32S3«335» 0.2031674<672307, 

4 0.23J49I5S6536J5, 0.2491 470458 1340, 0.0, 0 . 40484004765316E-1 , 

5 0.921 21499837 728C-1, 0.1 3887351 021 979, 0,1 78143980 761 95, 

6 0.20781604753699,0.22628318026290,0.23255155323087, 

7 0.351 19460331 7522-1,0.801 58087 1597608-1,0.12151857068790, 

8 0.15720316715819,0.18553839747794,0.20519846372130, 

9 0.21526385346316/ 

C 

IF  CLBAR.LC.01  SO  TO  400 
1FCLBAR.0T.14)  60  TO  400 
N 4 CLBAR*l)/2 
NEV  4 0 

IFCN.C0.LBAR/21  NEV  4 1 

00  300  J«lfN 

HBOC J)  • CHOCJ,LBAR) 

6XCJ)  • -CCA  C J , LBAR) 

IFCJ.EO.N.AND.NEV.ECI.O)  GO  TO  300 
M ■ LBAR  - J ♦ 1 
6X00  • -6XCJ) 

HGQCMJ  4 H60CJ) 

300  CONTINUE 
350  RETURN 
C 

400  MRITCC6, 1009)  LBAR 

1000  FORMAT  C29H0THE  VALUE  OF  LBAR  IS  INVALID/ 

1 8N0LBAR  4 14) 

LBAR  t 0 

RETURN 

ENO 
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• DECK  LIST  1 
SUBROUTINE  LIST1 

THIS  SUBROUTINE  PRINTS  AND/OR  CHECKS  CRITICAL  STRAINS#  STRESS  AND 
DISPLACEMENTS  for  The  multilayer  (NOfRVtt)  methoo, 

KZ»  PRINT  CODE 

0,  RETURN 

1#  COMPUTATIONS  ONLY 
2#  DON'T  CHECK  MAX  BUT  DO  PRINT 
3#  CHECK  MAX  AND  PRINT 
NOSE  • USE  CODE  FOR  SPATIAL  POINTS. 

0#  NO  USE. 

1,  PRINT  ONLY. 

2,  INTEGRATION  PURPOSES  ONLY. 

3,  PRINTOUT,  TOO. 

"CALL  CBLK1 

•CALL  CSLK6 
•CALL  CBLKR 
•CALL  CBLK10 

• CALL  CSLKM 
•CALL  CNOVA 
•CALL  CBLANK 
C 

OATA  FL1/2H  / , FL2/2H  •/ 

C 

IF  CKZ.EO.O)  GO  TO  1000 
IF(KZ.Eti.l)  GO  TO  10 
IF  INCALL.EO.O)  NRITE  (6,5000)  TIME 
NRITE  (6,5100) 

DO  5 I«l , MG 
M « HGM(I) 

DO  5 J*1 , MB 
N • NBN(J) 

IF  (MUSKJ,  I)  .EQ.O)  60  TO  5 

NRI TE (6, 5200)  M, N, UU( J , I ) , VV ( J, I) , MN( J , l ) 

5 CONTINUE 

IF  (NKP.EU.O)  GO  TU  1000 
IF(NPL'I.EQ.O)  NRI TE (6, 5600) 

IF (NPL  T.EO. 1 ) NRI TE (6, 5300) 
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00  750  1*1  * N6T 

X ■ X6(I) 

00  700  Jal.NBT 
NNU9E  • NU8E(J,I) 

IF  ( NNU9E.E6.0)  60  TO  700 

M B M ♦ 1 

IF  (NNU8E.LT, 2)  60  TO  700 

IF  (KZ.E0.2. AN0.NNU8E.EQ.2)  60  TO  700 

THSXB(J) 

EXX  S X 1 A ( M) 

ETT  a X2A (M) 

EXT  a XSA(M) 

XKXX  a X4A (M) 

XXTT  a X5A(M) 

XKXT  a XbA (M) 

00  650  Ral.NZP 
II  a NLZ(K) 

82  a ZCOO/A 

XI  8 EXX  ♦ 82*XKXX 
X2  a ETT  ♦ 82*XKTT 
XI  a EXT  ♦ 82*XKXT 

816XX  a BXL(IZ)* (XI  ♦ THNU(II)*X2) 

SI6TT  8 BTL ( 1 1 ) * (Xi  ♦ XXNU(1I)*X1) 

SI6XT  a 6XT (1I)*XJ 
C PRINCIPAL  BTRE99E9. 

816  a SORT ( ,25* (8I6XX  • SI6TU**2  ♦ SI6XT**2) 

81  a (&I6XX  ♦ 8I6TT) *.5 

8161  a 81  ♦ 816 

8162  a 91  • 916 
600  FLA 6 a FL1 

IF  (8I61.6T.9AT(II))  FLA6  a F(.2 
IF  (8I62.LT,«9AC(II))  FLAG  a Fl2 

NR1T8  (b* 5400 ) X,TH,ZC(K),X1, X2, XI, 8I6XX, 9I6TT, SIGXT.FLAG 
b50  CONTINUE 
700  CONTINUE 
750  CONTINUE 

IF  (KZ.LT.2)  60  TO  1000 
C 

NRITE  (b# 6000) 
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DO  7/0  l * l * MB AH 
X « X5(I) 

00  770  Jsl'NBAH 
NBC  ■ NBU3E(J,I) 

IP  (NBC.EO.O.OR.NBC.GT.lOO)  BO  TO  770 

Y * X8( JJ 

IF  (NBC.GT.Q)  BO  TO  760 
NC  = -NBC 

WRITE  (6,6100)  X,Y,VRX(NC),ENX(NC) 

BO  TO  770 

760  WRITE  (6,6100)  X , Y , VHT (NBC) , ENT (NBC) 
770  CONTINUE 

IF  (NBND.EB.l)  60  TO  790 
WRITE  (6,6200) 

00  780  1*1, NRC 
X * XQ  ( 1 ) 

IF  (1.GT.2)  X * XG (MB AH ) 

Y « XBU) 

IF  (UEQ.2.0R.I.EQ.8)  Y • XB(NBAR) 
WRITE  (6,6100)  X,Y,RH(I) 

780  CONTINUE 

790  IF  (NPlT.EO.O)  WRITE  (6,5700) 

IF  (NPLT.EQ.l)  WRITE  (6,5500) 

KKK  * 0 

00  900  1*1, NGT 
X * XBU) 

00  800  J* 1 , NBT 
NNUSE  « NUSE ( J, I) 

IF  (NNUSE. EO.O)  SU  TO  800 
KKK  * KKK  ♦ 1 
IF  (NNUSE. LE. 2)  60  TO  800 
IF  (NU.EU.O)  PPP  • P (KKK) 

UF  s A*U (KKK) 

VF  * A* V (KKK ) 

MF  * A*W(KKK) 

WRITE  (6,5600)  X , XB ( J) , UF , VF, WF, PPP 

800  CONTINUE 
900  CONTINUE 
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1000  return 

c 

5000  FORMAT  (//IH0,JNT  a,C15.7,4H  9EC) 

SI 00  FORMAT  (1M0,2X,5MBAMMA, 3X, 4MBCTA, BX, 3HUR8, 1 JX, 3NVR9, 1 3X , 3HNR8) 
5200  FORMAT  (IS, XT, 2X, 111  A. 7) 

5300  FORMAT  (1H0,3X, 1MX,JX,5H  BETA, 5X, IHZ, 10X, 3HEKX, 12X, 3HETT, 12X, 

1 3HEM T , 4X  » 0H816MA  XX, 7X, 8H91GMA  TT , 7 X , 8H9I6MA  XT/ 

2 iX,4M(IN),2X,5H(0t8),3X,4M(lN),6X,7M(IN/IN),9X,7H(IM/IN), 

3 BX,?M(IN/IN),4X,§H(P8I),10X,5H(P8I),10X,5H(PSI)) 

5400  FORMAT  (lX,F9t2,2X,F8.2, lX,F7,4,6611,6r2X,A2) 

5500  FORMAT  (1H0,7X,5HX(IN) ,9X,  10H  BETA (DC6) » 8X, SHU ( IN) , 

1  l 1X,5HV( IN) , t 1 X, SHM ( IN) , OX, 14HPRE99URE  CPSI)) 

5600  FORMAT  (1X,6E16.7) 

5700  FORMAT  ( 1M0, 7X , 5HX (IN) , 12X, 5MT ( IN) , 1 OX, 5HU( I N)  , 

1 1 IX, SHV ( IN) , 1 1 X, SHM ( IN) , 6X, 14HPRE88URE  (P8I)) 

5600  FORMAT  ( |HB, 3X, |HX,6X, 1HY , 6X, IHZ, )0X , 3HEXX, 12X , 3HETT , 12X , 

1 3HEXT , 4X , 0H8I8MA  XX, 7X, 0H8I6MA  TT, 7 X, 8H9I6MA  XT/ 

2 2X,4M(|N),1X,4H(IN),3X,4H(IN),6X,7H(IN/IN),BX,7HUN/IN), 

3 eX,7N(IN/IN),9X,ftH(P8I),10X,5H(PtI),10X,5H(P8I)) 

5900  FORMAT (27M00CFLECT10N  AT  CENTER, IN  a E15.8) 

6000  FORMAT  (51H0REACTI VC  F0RCC9  PER  UNIT  LENGTH  ALONG  E06E  (L8/IN)/ 
1 6X,|HX,tX,lHT,9X,lHV,14X,lHN) 

6100  FORMAT  (2F«,),2E15.6) 

6200  FORMAT  (33H0REACT1 VE  FORCES  AT  CORNERS  (LB9)) 

ENO 
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t nno onooonoonno 


• OECK  LIST? 

8UBR0UUNE  LIST? 

THIS  SUBROUTINE  PRINTS  AND/OR  CHECKS  STRAINS#  STRESSES#  AND 
DISPLACEMENTS  POM  THE  SINGLE  LAYER  (N0CRW2)  METHOD* 


KZ-  PRINT  CODE 
0#  RETURN 

1#  COMPUTATIONS  ONLY 
?#  DON'T  CHECK  MAX  BUT  DO  PRINT 
3#  CHECK  MAX  ANO  PRINT 
NUSE  • USE  COOE  POR  SPATIAL  POINTS. 
0#  NO  USE* 

I#  PRINT  ONLY. 

?#  INTEGRATION  PURPOSES  ONLY. 

3#  PRINTOUT#  TOO. 


CALL  CBLK1 
•CALL  CBLK7 
•CALL  CBLKB 
•CALL  C6LK9 

• CALL  CBLMO 

• CALL  CBIMA 
•CALL  CNOVA 
•CALL  CBLANK 
C 

DATA  PL1/2H  / # PL?/?H  •/ 

IP (KZ.EQ.O)  GO  TO  1000 
IP(KZ.EG.I)  GO  TO  10 
IP  (NCALL.EO.O)  NRITE  (6#  5000)  TIME 
NRITE  (6,5100) 

DO  5 I •! # MG 
M • M6M(I) 

00  5 J*1 # MB 

IP  (MUSE(J#  D.EU.O)  GO  TO  5 
N • NBN(J) 

*RITE<6#5?00)  M#N#UU(J#I)#VV(J#I)#NW(J,I) 
5 CONTINUE 

IP  (NKP.EO.O)  GO  TO  1000 
IP(NPLT.EO.O)  NRITE(6#5B00) 

IP(NPtT.EG.l)  NRITE (6# 5300) 
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10  Mao 

00  300  lal'NST 
X a *•<!) 

00  200  jal,NBT 
MNU8E  a MU6E ( J, I ) 

IF  (MMU3E.EO.0)  80  TO  200 
M a M ♦ 1 

IF  (MMU3E.LT, i)  60  TO  200 
IF  (RZ.EB.i. AND.NNU8E.E9.2)  60  TO  200 
JI  a L0AR* (M-l ) 

THBXI(J) 

Exx  a X 1 A ( M) 

ETT  a XtA(M) 

EXT  B X 3A ( M) 

XKXX  B X4A (M) 

XKTT  a X9A(M) 

XKXT  a XOA (M) 

DO  100  KRal, 2 
X a IMZ(RR) 

L « JI  ♦ R 
32  a ZA(RR) 

31  a 32 

IF  (MCLP.EQ.2)  81  a ZF(R) 

flamfli 

IF  (MfLF.EO. 1 . ANO.EFBO(L) « 8 T, 81602)  FLA6  > FL2 

XlafRX  ♦ SZftXRXX 

X2a  ETT  ♦BiaKRTT 

X3a  IRT  ♦ S2*XRXT 

IF  (MELF.E8.2)  60  TO  15 

8181  a 81 A (M)  ♦ 32*84A(M) 

3162  • 82A (M)  ♦ 82*85A(M) 

3183  a 33A (M)  ♦ 32«36A(M) 

60  TO  00 

15  8161  • 3XX (L) 

8182  a 8TKU 
8163  • 8XT (L) 

00  IF  (MMU8E.E6.2)  60  TO  100 
IF  (R2.E0.1)  60  TO  100 

MRITC  (0,5400)  X,TH,*B(RR), XI, X2,X3, 8161, 3162, 3183, M(U,FLA6 
100  COMTIMUE 


LIST2 


166 


aoo  CONTINUE 
300  CONTINUE 

IF(KZ.LT.i)  60  TO  1000 


WRITE  (6,6000) 

00  370  I*1,MBAR 
X * X6(I) 

00  370  J*l , NBAR 
NBC  ■ NBUSE( J, I) 

IP  (NBC.EQ,0.OR,NBC.6T,100)  60  TO  370 

Y * XB(J) 

IP  (NBC.6T.0)  60  TO  360 
NC  • «NBC 

WRITE  (6,6100)  X#Y,VRX(NC),ENX(NC) 

60  TO  370 

360  WRITE  (6,6100)  X, Y, VRT  (NBC)  , ENT  (NBC) 
370  CONTINUE 

IF  (NBNO.EO.l)  60  TO  390 
WRITE  (6,6200) 

00  300  1*1 , NRC 
X a X6(l) 

IP  (I.6T.2)  X ■ XG(NBAR) 

Y • X B ( 1 ) 

IP  (I.EQ.2.0R.I.EQ.4)  V a XB(NBAR) 
WRITE  (6,6100)  X,Y,RR(I) 

380  CONTINUE 

390  IP(NRUT.EU.O)  WR1TE(6, 5700) 
IP(NRUIT.EO.I)  WRITE(6,5500) 

KKKaO 

DO  500  1*1 , N6T 
X a XB(I) 

00  400  J*1 , NBT 
NNU8E  • NU8E( J, I) 

IP  (NNUSECE6.0)  60  TO  400 
KKX  • KKK  ♦ 1 
IP  (NNUSE.LE.2)  60  TO  400 
IP  (NU.EQ.O)  PPP  a P (KKK) 

UP  * A*U(KKK) 

VP  a AaV(KKK) 
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«F  • AftM(KKK) 

NRITE-  (6,5600)  X,X»(J) ,UF, VF,RF,PPP 
400  CONTINUE 
500  CONTINUE 
C 

1000  RETURN 
C 

5000  FORMAT  (//1M0,3HT  *,E15,7,4N  SEC) 

5100  FORMAT  (1M0,2X,5HCAMMA,3X,4HBETA,8X,3HURS, 13X,3HVRS, 13X,3H*RS) 
5200  FORMAT  ( IS, 17, IX, 3116. 7) 

5300  FORMAT  (1MQ,3X, 1HX,3X, 5H  SET A, 5X , 1 HE, 1 OX , 3HEXX , 12X , 3H|T T , 12X, 

1 1HEKT, OX, SH310MA  XX, 7 X , 8H3I8MA  TT, 7X, 8HSI6MA  XT , 3X, 6HREGIQN/ 

2 2X,4H(IN),2X,5M(0E0),3X,4H(IN),6X,7M(IN/lN),«X,7H(IN/IN), 

3 SX,7H(IN/IN),0X,5H(P8I),10X,5H(P8l),10X,5H(PSI)) 

5400  FORMAT  ( 1 X , F5. 2, 2X , F5. 2, 1 X , F7 ,4 , 6E1 5,6, 14, A2) 

5500  FORMAT  ( 1H0,  7X,  5MX  (IN) , OX,  10H  SET  A (DEC) , OX,  5HUUN)  , 

1 1 1 X , 5MV (IN),11X,5HM(IN),6X,1 4HPRE38UKE  (P91)) 

5600  FORMAT  (1X,6E)6.7) 

5700  FORMAT  ( 1M0 , 7X,5MX (IN) , 1 2X , 5MT ( IN) , 1 OX , 5HU ( IN) , 

1 1 IX, SMV ( IN) , 1 IX, 5HM( IN) , 6X, 14HPRE88URE  (PSI)) 

5800  FORMAT  (1H0, 3X , 1HX ,6X , 1HT , 6X , 1HZ , 1 OX, 3HEXX, 12X , 3HET T, 1 2X , 

1 3hEXT,0X,6M3ISMA  XX , 7 X , SRSIGMA  T T , 7X , 6H8IGMA  XT, 3X , 6HKEGI0N/ 

2 IX, 4H( 1 N), SX, 4M ( IN), 3X, 4H( IN), 6X, 7H( IN/ IN), «X,7M( IN/ IN), 

3 SX,7H(IN/IN),9X,5H(P81),10X,5H(PSI),10X,5H(P8I)) 

6000  FORMAT  (51HQREACT 1VE  FORCES  PER  UNIT  LENGTH  ALONG  EDGE  (LS/IN;/ 
1 6X, 1 HX , 8X , 1 HY , OX , 1 HV, 1 4X , 1 HN) 

6)00  FORMAT  12F0. 3, 2115,6) 

6200  FORMAT  (33H0RCACTIVE  FORCES  AT  CORNERS  (LBS)) 

END 
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• DEC*  PINIT 

SUBROUTINE  PINIT(H) 

•CALL  CNOVA 
•CALL  CLOAO 
•CALL  CBLftl 
•CALL  CBLK6 
•CALL  CBLMNA 

DIMENSION  001(3), 002(3) 

DATA  001/87, E-6,  90. E-6,  84.E*6/,  002/0.0*86,  0.1127,  0.1275/ 

C 

IF(M.|0.1)  60  TO  200 
IF (KD8.EQ.2)  60  TO  150 

STATIC 

RE AO (5, 2000)  PS 
*RITE(6,2200)  PS 
NU«1 
PPPaPS 
150  RETURN 

OYNAMIC 

200  IF (KDS.EO. 1 ) 6U  TO  400 
RE AO  (5,1000)  NLOAO 
NRITE  (6, 2«00)  NLOAO 
60  TO  (600,800,250,500),  NLOAO 
250  READ (5, 2000)  PPl , PPO, TTQ, TPR1ME, AA, ANN 
NRITEC6, 2300)  PPl , PPO, TTO,TPRIME,AA,ANN 
NU«1 

IF (T PR  I ME, EQ .0.0)  60  TO  300 
PPRIM€»PPO*(1,0  • TPRIME/YTO)**ANN 
PPRIME  • PPRIME*EXP(*AA*YPRIME/TTO) 
TTl«TPRlMt*PPl/(PPl«PPRXME) 

OTT 1*1 ,0/TTl 
300  OTTOai.O/TTO 
AZaA AaOTTO 
400  RETURN 
C 

500  READ  (5,1000)  NTIME 
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mo  (5*2100)  (TT(1),PT(I),I»1,NT1ME) 

MRITE  (6*1500)  NTIME,  (TT(I),PT(I),I»1, NTIME) 

NU  • 1 
JL  • 2 
RETURN 

C (RUIN  SYMMETRIC*  NONUNIFORM  LOAO  ON  FLAT  FLATS. 

600  VS  • 5.SSEA 

READ  (5*2000)  ZEE* FHI 
MRITE  (6*2600)  ZEE, PHI 
XO  ■ XO(NS) 

YO  ■ XI (NS) 

00  700  I*1*NST 
x ■ xs(i)  • xo 
00  700  J*l  * NBT 
V * XS(J>  • YO 

700  RA ( J* I ) • SORT (X**2  t ♦ ZEE*»2) 

IFHI  • 1 

IF  (FHI.0T.O.O)  IFHI  • 2 
IF  (PHI. ST, 10.0)  IFHI  « 3 
U1  ■ BQl(IPHI) 

02  ■ a02(IFHI) 

NU  • 0 
RETURN 
C 

600  READ  (5*1000)  NPX*NPY*NT1ME 
REAO  (5*2000)  DTIM 
MRITE  (6*2700)  NPK* NPY* NTIME* OTIM 
REAO  0*2000)  (XP(I)*J«1* NPX) 

NRXTE  (6*1100)  (XP(X ) * 1*1  * NPX) 

REAO  (5*2000)  (YP(J)  , JM,NPY) 

MRITE  (6*1200)  (YP(J)*J«1*NPV) 

MRITE  (6*3300) 

00  SIO  I«1*NPX 

REAO  (5*2000)  (0ET(J,I),J«1,NPY) 

020  MRITE  (6* 2S00)  (OET( J, I) , J»1 , NPY) 

MRITE'  (6*2000) 

00  S40  111, NPX 
00  S40  J«1,NPY 

REAO  (5*2000)  (PRT (X, J, I) * X«1 , NTIME) 

8*0  MRITE  (6*1000)  (PRT(K, J, I) * X«1 , NTIME) 
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C SPATIAL  INTERPULATION»EXTRAPOLATION.  INDICES  ARE  LONER  BOUND* 
DO  900  1*1, NOT 
DO  860  III  • 1,NPX 
IF  (XPUin.CT.XSU))  DO  TO  880 
860  CONTINUE 
III  * NPX 

880  IF  (I1I.GT.1)  HI  • HI  • 1 

0X1(1)  a (Xfi(I)  • XP(UI))/(XP(HItl)  • XP(III)) 

900  1X1(1)  a HI 

00  960  J a 1 , NBT 
00  920  JJJ  • l » NPY 
IF  (YP(JJJ).6T.XB(J))  60  TO  940 
920  CONTINUE 
JJJ  » NPT 

940  IF  (JJJ.6T.1)  JJJ  a JJJ  - 1 

DTl(J)  a (XBC J)  • YP(JJJ))/(YP(JJJM)  - TP(JJJ)) 

960  JVJ(J)  a JJJ 
NU  a 0 

00  980  I*lfNPX 
DO  980  Jal'NPY 
980  JLTIJrl)  • 2 
RETURN 
C 

1000  FORMAT  (6112) 

2000  FORMAT (6F 12* 1 ) 

2100  FORMAT  (2F12.1) 

2200  FORMAT (2RH08TATIC  PRESSURE*  PSI  a E1S.6) 

2300  FORMAT (23H00YNAMIC  LOAD  CONSTANTS/ 


11H 

PP1 

a 

£15*6/ 

11H 

PPO 

a 

115.6/ 

UN 

TTO 

a 

€15.6/ 

1 1 H 

TPRIME 

a 

£15.6/ 

UN 

AA 

a 

E15.6/ 

1 1 H 

ANN 

a 

€15.6) 

2400  FORMAT  (21H0DVNAMIC  LOAD  OPTION  14) 

2500  FORMAT  (18H0NUMBER  OP  TIMES  • I4/2BH  TIME*  SEC  PRESSURE*  PSI / 
1 (2E1  5*6) ) 

2600  FORMAT  (41H00VNAMIC  LOAD  CONSTANTS  • FLAT  PLATE  ONLY/ 

1 14H  ZEE  (IN)  a E15.6/  14M  PHI  (DEC)  * £13,6) 

2700  FORMAT  (23H0DYNAMIC  LOAD  CONSTANTS/ 
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1 1IH  HP*  ■ I5/12H  MPY  • IS/ 

2 1 iH  NTIME  • IJ/12H  DTIM  ■ E15.6) 

2*00  FORMAT  (9K,9El9,t) 

2900  FORMAT  (1 2H0PRC8SURES  ») 

1000  FORMAT  (9X,*£l5tS) 

3100  FORMAT  (19HOX-POSITION8  UN)  »/  (5X , 5E 1 5. 0)  ) 

3200  FORMAT  (26HOY-POSI TI0N8  (IN  OR  DEG)  ■/ (5X f SE 19,6) ) 
3300  FORMAT  (20H00ELA Y TIMES  (SEC)  «) 

EMO 
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• DEC*  PRESS 

subroutine  press 
•call  cnova 
•CALL  CLOAO 
•CALL  CBLK1 
•CALL  CBLKS 

• CALL  CBLANK 
C 

IF (NCALL.GT.O)  SO  TO  1000 
ZZ«  l.O/RTRIALU) 

SO  TO  (400,600,50,220),  NLOAO 
50  IF ( TIME, GEa  7 PRIME)  SO  TO  100 
PPP4ZZ«PP1*(1,0  • TIME*0TT1) 

IF(PPP.LT.O.O)  PPP»0,0 
60  TO  1000 

100  IF(TIM£,SE,TT0)  60  TO  200 

PPP4PP0*(1,0  • TlME»OTT0)**ANN 
PPP4ZZ*PPP*EXP(-AZ*TIM£) 

60  TO  1000 
200  PPPlO.O 

60  TO  1000 

c 

220  00  240  J»JL, NTIM6 

IF  (TIMi.LE.TTCJ))  SO  TO  2S0 
240  CONTINUE 

WRITE  (4,250)  TIME, TT(NTIME) 

250  FORMAT  ( J2H  *•  WARNING  « TIME  EXCEEDS  TABLE,  2E15.6) 
J • NT  I ME 
2bO  JL  « J 

PPP  « PT(J.l)  ♦ (TIME  • TT(J.l))*(PT(J)  • PT(J-l))/ 

1 (TT(J)  - TT (J»l) ) 

PPP  « ZZ*PPP 
60  TO  1000 
C 

400  OUM  4 TIME  ♦ ZEE/VS 
TO  • 81  * Q2*DUM 
PM  ■ 464,fS*0UM**(»,29) 

4 4 0 

00  600  141, N6T 
DO  600  Jsl,NBT 
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IF  (NUSE(J,I).E8,0)  60  TO  600 
K » K ♦ t 
RRR  ■ RA(J,I) 

T8TAR  • (RRR  - ZEE)/V8 
PPP  • 0,0 

IF  (T8TAR.8T.TIME  ♦ l.OE-12)  SO  TO  500 

IF  (TIMf*T»TAR.ST.TO)  00  TO  500 

FAC  • (TITAR  • TIME) /TO 

CALPH  • ZEi/RRR 

PPP  • PM*CALPH»(1.0  ♦ FAC) 

500  P(K)  • PPP*ZZ 
600  CONTINUE 
80  TO  1800 
C 

C INTERPOLATE  ON  TIME, 

800  00  6b0  I ■ l » HP X 
00  860  JMiMPY 
PPP  • 0.0 
OETT  • OET ( J»  I ) 

IF  (TIME.LT.OETT)  60  TO  860 
JL  • JLT(J,1) 

00  818  K«JL*NTIME 
KK  = K 

IF  (TIMC.LC.DETT  ♦ DTIM*FLQAT (K-l ) ) 60  TO  840 
820  CONTINUE 

JLT l J» I)  • NTIME 
PPP  • PRT (NT 1ME. J. I) 

60  TO  860 
840  JL  * KK 

PI  • PRT(JL*1,J,I) 

U • OITT  ♦ 0TIM*FL0AT(JL«2) 

PPP  » PI  ♦ (TIME  • T 1 ) * (PR T(JLsJ*I)  • PD/DTIM 
JLT (Jp I)  m JL 
860  PRTT ( J, I)  • PPP 
C INTERPOLATE  SPATIALLY. 

K ■ 0 

DO  880  1*1, N6T 
II  • 1X1(1) 

OX  • 0X1(1) 

DO  880  Jal'NBT 
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IF  (NOSE ( J * I ) , EQ, 0)  60  TO  ISO 
K • K ♦ 1 
JJ  • JTJCJ) 

0»  • OU(J) 

PI  » PRTT(JJ,II)  ♦ DY*  (PRTT(JJ*1 # II)  . PRTT(JJ,IJ)) 

Pi  ■ PRTT(JJ,II*l)  ♦ 0T»(PRTT(JJ*1,  IIM)  • PftTT(JJ,IIM)) 
PPP  • PI  ♦ OX* (P2  • PI) 

PCX)  s PPP*ZI 
8ttO  CONTINUE 

1000  RETURN 
END 
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ftOCCK  PROP 

SUBROUTINE  PROP 

C 

•CALL  Ft«»T 
• CALL  CSLK  t 
•CALL  COLKi 
•CALL  CBLK1 
•CALL  CSLKA 
•CALL  CBLRS 
•CALL  CBLK7 
•CALL  CBIKS 
•CALL  cblko 
•CALL  CBLKIO 
•CALL  CBLK1I 
•CALL  CBLA1J 
•CALL  CBIK14 
•call  CNOOA 
•call  cblan* 
c 

IP(NCALL.EO.O)  60  TO  5700 
IP(NCAUtCQtl)  60  TO  5450 
PI>1»  1 41 99205IS00A 
CRRIsf  » 0E»5 
IP1R8TB0 
JFIHST  • 0 
NO  • 1 
C 

CALL1  OBCTt 

C 

call  obetb 

1F(K|RR.6T.0)  60  TO  6300 

C 

CALL  OBETS 

C 

00  5300  M«1,M6 
00  5300  NM'MB 
U1(N,NOM.O 
Vl<N,N)»0.0 
* 1 IN,  *0*0.0 
5300  CONTINUE 
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00  5000  J*1,NY2 
XX(J)«0,0 
5400  ERRU)«ERR6 
60  TO  5900 
C 
C 

5450  IF (K0S,£0,i)  60  TO  5900 

IF  (PPP.EB.O.O)  60  TO  5900 

STATIC  SOLUTION 

NNELP  • NELF 
NELP  » 1 
KZ  i 2 
NCQUNT  s 0 
NOKKO 
NTRaO 

MTRMAX  • 10 
CPTIM1  • 8EC(0UM) 

5500  CALL  DERV2 

CALL  RELAXP(NY2,YY,XX,ERR,N0X,NDBUG, NGOUNT) 

ICOUNT  • ICOUNT  ♦ 1 

NTRaNTR+1 

MTR»(NTR-2)/lNY2*l) 

IF  IMTR,E0.(NTR*J)/(NY24-1))  60  TO  5550 
IF  (NDBU6CE6«0)  60  TO  5500 
IF  INDERV.EQ.l)  CALL  LXST1 
IF  (NOERV.EO.2)  CALL'  LI8T2 
5550  IF (MTR,6T,MTRMAX)  60  TO  6200 
IF(NOft.EQ.Q)  60  TO  5500 
IF(N0K*6E.2)  60  TO  6300 
MTRa(NTR*l)/(NY2M) 

NKLP  • NNELP 
J»1 

00  5600  Mal,M6 
00  5600  N«1,M6 

IF(MUSE(N,M) •E6a0)  60  TO  5600 
U1 (N,  4) aXX  ( J) 

V1(N,NI«XX(J«M6M6) 

Nl(N,N)aXX(J«H6MB2) 
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J»Jtl 

5600  CONTINUE 

IF (KOS«f B«1  . AND*  NOSUO.EQ.O)  SO  TO  5900 
CPTJM5  a fECIOUMJ 
CPTaCPTIMJ*CPTIMi 
nRITI16«SI50)  MTR*  CPT 
IF  (NOCRV,E0.U  CALC  LI8T1 
IF  (NOCRV tE0«2)  CALL  IIST2 
00  TO  5900 
C 

C DYNAMIC'  RESPONSE 
C 

5700  iFIRSTao 

CPTIMI  • SEC(OUH) 

Jal 

KC  • 0 

DO  5750  Mil, MO 
00  57S0  Nat i MS 

IF IMUSEINf M) tE8.0)  60  TO  5750 
XX(J)aUl(N,M) 

XX( J9M6MS)  a VI (Nf M) 

XRU»M6MS()aNl(NtM) 

VXOUIM, 

VXO(J«MOMB)aO. 

VXOU»MSMS2)ao, 

JaJ»| 

5750  CONTINUE 
SMAXaO.O 

OPRTa-0.5*OELTlM 

TIMEaO.O 

KHIH  9 1 

TF  a TSTOP  ♦ v9*0CLtTIM 
5760  CALti  OCRVe 

CALL1  HIM  (KHIM'NVifDELTIM'TlMlfVXO'XXrVV) 

ICOUNT  a 1C0UNT  ♦ t 

IF  (X|ftR.BT«0)  60  TO  6400 

IF  (TIME«LTVTF)  60  TO  5760 

MRITC<  (6»  1 1400)  TIME 

IFCKOAM.CO.f)  60  TO  5S00 

CRIT (l)aSMAX 


o r»  o n n n n 


WRi  ft (b, 8400)  NTRIAL,NCASE,RTRIAL(1),SMAX,TMAX 
IF(NELP  .68,4)  WRITE (6, 8900)  NRE6 
IFlNTECU.EO.l)  WRITE<6,6500) 

IF (NTEC0.E0.2)  WRITE (8.6600) 

5800  CPTIM3  « SEC (OUM) 

CPT  » CPTIM3  * CPTIMI 
WR I TE(6» 6000)  CPT 
IF ( JFSRST.EQ.O)  60  TO  5*00 
N8«0 

00  5650  L»1,LMAX 
IF(KY(U,6T,1)  N8»NS*1 
<3850  CONTINUE 

*RITS(6,*200)  N8,LMAX 
5900  RETURN 

ERROR  MES8A6E8 

6200  WRITE (6, 6*50)  MTR 
6300  KERR  » 1 
TIME  • 0, 

6400  WRITE  (6,11500)  TIME 
RETURN 

FORMAT  STATEMENTS 

6350  FORMAT (26H1RE6ULT8  OF  STATIC  PRELO AO /2 OHO NUMBER  OF  TRIALS  • 15/ 

1 ISM  NET  CP  TIME  ■ P11.S) 

6400  FORMAT (1TM0RE6ULT8  OF  TRIAL1S,6H  OF  CA8EI3/21H  RANSE*  FT  ■ 

1 E15.6/21H  CRIT(l)  • E15.6/I1H  TIME,  SEC  ■ US. 6 

2) 

6650  FORMAT  (I5H0T00  MANY  TRIALS  IN  STATIC  SOLUTION/TH  MTR  • 14) 

6500  FORMAT (10M  TENSION) 

6600  FORMAT (14H  COMPRESSION) 

6*00  FORMAT (25H  ELASTIC-PLASTIC  RE6I0N14) 

9000  FORMAT  (33M0NET  CP  TIME  FOR  RESPONSE,  SEC  • F10.3) 

*200  FORMAT  (1M0,I6,3H  0F,I4,15H  POINTS  VIELOEO) 

11400  FORMAT (1H6//42H  NORMAL  OEPROP  STOP  CONDITION  AT  T,  SEC  » E14.6) 
11500  FORMAT  (IN  //31H0DEPR0P  IS  ABORTED  AT  T,  SEC  « 114. 6) 

END 
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•DICK  REIT 

SUBROUTINE  RE1T(I* J.K) 

C CONFUTE  REACTIVE  FORCES 

•CALL  CSLM 

•CALL  COINS 

•CALL  CSLN7 

•CALL  CSLR10 

• CALL  CBLNU 

•CALL  CBLKtA 

•CALL  CBLANN 

NBC  a NBUOE(J,I) 

IF  (NOERV. E0.2)  60  TO  BOO 
C NOERV  a 1, 

IF  (NBC. 6T. 100)  60  TO  100 
IF  (N0C.6T * 0)  60  TO  100 
NBC  a *NBC 

VRX(NSC)  a C1M666(NBC)  ♦ C2MBBB  (NBC) 
ENK (NBC)  a CM11*EXX  ♦ CM12*ETT 
IF  (I.EO. 1 ) VRX(NBC)  a •VRX(NBC) 

60  TO  1100 

J 

100  VRT(NBC)  a CSaMBBS(NBC)  ♦ CR*A66B(NBC) 
ENT (NBC)  a CM22AETT  ♦ C*12*E«X 
IF  U.E6.1)  VRT(NBC)  • -VRT(NBC) 

60  TO  1200 
C 

200  NBC  a NBC  • 100 

RR (NBC)  a C5«nBB(K) • (•! a0) •• (NBC/t) 

60  TO  12B0 

NOERV  a 2, 

ELAITIC-FLAOTIC, 

BOO  BUM  a 0.0 
SUM1  a 0. 

JI  a lBAR*(K»1) 

IF  (NBC. 8T. 100)  60  TO  1000 
IF  (N6C.6T.0)  BO  TO  800 
NBC  a »NBC 
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Jll  = LBAK»(K*2) 

<3  • I*NBT  ♦ J 

IF  U.6T.1)  Kl  9 (1*2) *NBT  9 J 
Jll  « LBAR* (K2  • l) 

QELXX  • DELX 

IF  II.EO.I)  DELXX  • *DEIX 
DO  700  KK  s 1,LBAR 
L 9 JI  9 KK 
LI  9 Jll  ♦ KK 
Li  • J12  ♦ KK 
SI  9 SXXIL) 

611  • SXX(Li) 

S3  9 SXTCU 

631  9 SXT(Ll) 

SUM  9 SUM  ♦ ((611  • 6D/DELKX  ♦ 2.0»(631  • G3)/DELT)» 
1 6X (KK) *NG0 (KK) 

700  BUM1  9 6UM1  9 H60(KK)*61 
VKX(MBC)  • C5*BUM/i,0 
IF  (1,10,1)  VOX (NBC)  9 •VRX(MBC) 

ENX(NBC)  9 H9BUM1/2, 

GO  TO  1200 
C 

800  Kl  9 K • 2 

IF  ( J,EO, 1 ) Kl  9 K 

Jll  9 lbar*ki 

DELT  T 9 OILT 

IF  (J.EQ.l)  OELTT  9 •OCLT 
Jit  « LBAR*((1«2)*NBT  ♦ J*|) 

DO  900  KK  9 1* LBAR 
L i JM  KK 
Li  9 Jll  ♦ KK 
Li  9 J 12  9 KK 

62  9 BTT(L) 

622  9 STT(Ll) 

63  9 8XT (L) 

632  9 8XT (L2) 

SUM  9 BUM  9 ((622  • 62) /DELiTT  9 2,0*(632  • 63)/DELX)* 
1 6X(KK)»M60(KK) 

900  6UM1  9 8UM1  9 H60 (KK ) *62 
VRT(MBC)  9 C»*»UM/2,0 
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IP  IJ.tU.l)  VRT(NSC)  • •VRT(NBC) 
(NT (NSC ' • HftSUMl/2,0 
60  TO  1200 
C 

1000  NSC  ■ NSC  • 100 

00  1100  KKai'LSM 

1 • JI  ♦ KK 

1100  SUN  • SUN  ♦ 6X(KK)*NS0(KK)«SXT<1) 
SR (NSC)  0 CfaSUN* (*1,0) *• (NSC/I) 
1200  Nf TURN 
(NO 
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• DEC*  «ELAXP 

SUBROUTINE  RELAXP  (NEQ, RES, X , ERR, NOR, NPRINT , NCQUNT) 

DIMENSION  RES (l)fX(l)f ERR ( 1 ) 

• CALL  cblku 

DATA  CON/SOOO,/ 

IF  (NPRINT. NE. 2)  GO  TO  40 
RR1TE  (6,50)  (X(N),RES(N),N*l,NEa) 

50  FORMAT  (1H  ,//«X,6H  TRIAL  X, 1 OX , 7HRCS1DUI/ (4X , (4 J,6, 4X,E1S,6) ) 
40  IF  (NCOUNT.EB.O  ) BO  TO  10 
IF  (NCOuNT.lE.NEO)  80  TO  14 
NCOUNT*0 
N0*0 
NN8*0 

00  2 1*1, NEO 

OX  « X(I)  • PX(I) 

IF  (A88(0X).IE,ERR(I))  NO  s NO  ♦ l 
IF  (ASI(RE8(I))*ABS(PRE8(I)),6T,0,0)  NNG*NNQ»I 
2 CONTINUE 

IF  (NO.EQ.NEQ)  GO  TO  100 
IF  (NNB.E8.NEB)  GO  TO  101 
10  00  4 1*1, NEO 

XX1(I)«X(I) 

RRE8(I)*RE8(I) 

0£LX(I)*AB8(0.0001*X(I)) 

IF  (OELX(I)vLT,ERR(I>)  DELX(I)*ERR(1) 

4 CONllMUt 
GO  TO  6 

14  00  7 MM»1 , NEQ 

7 XRES(MM,NCOUNT)  » (RES (MM)  • RREB(MM) ) / DELX (NCOUNT ) 
X(NCOUNT)bXXKNCOUNT) 

IF  (NCOUNT, EQ.NEQ)  60  TO  S 
6 NCOUNTBNCOUNWl 

X (NCOUNT) «X (NCOUNT )t0ELX (NCOUNT) 

NOK«0 
RETURN 

6 00  20  laifNEQ 
20  8I6X(I)*»RRE8(1) 

CALL  SOLVE  (XREB, NEB, 147,0, (F,DET,S16X) 

IF  (NC9.EQ.0)  80  TO  15 
PROP*! ,0 
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00  II  lai'NCO 

IP  ( A00(8]CX  d) ) •LTcCON*ENR(  I ) ) 00  TO  IS 
XMOf  • COM*fftR (I)  / AMOIBX  ( I) ) 

IP  (XMQPck.T«MOP)  PROPbXMO? 

II  COMTZMOC 

00  t»  UtfMtO 
KU>««Rld)+tlSXd>*PROP 
PXd)MMd) 
l i PRCO(X)MRCOd) 

NOKBO 

NCOUNTbNEOpI 
RtfURN 

100  NQIUI 
00  TO  d 

101  RRITI  Ut99) 

99  PORNAT  (UNO 
19  NOR  • i 
11  NCOMNUO 

RETURN 
(MO 


SOLUTION  OIVERSIMC  IN  RELAX?) 


♦DEC*  SEC 

FUNCTION  SEC  (DUM) 

C FIND  ELAPSED  CF  TIME. 

C 

CALL  SECOND  (SEC) 

SETUSN 

END 
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• DICK  sigma 

SUBROUTINE  SIGMA  (I,J,M) 

C 

C THIS  s/r  determines  the  stress-strain  RELATIONSHIPS  pqr 
C ELASTIC  AND/OR  PLASTIC  RESPONSE* 

C SUBROUTINE  COMPLETELY  REVISED  MARCH*  197b. 

C 


c 

A 

• INOEX 

Of 

THE 

INTEGRATION 

POINT 

IN 

THE 

Z DIRECTION. 

c 

1 

- INOEX 

Of 

THE 

INTEGRATION 

POINT 

IN 

THE 

BETA  DIRECTION. 

c 

J 

• INOEX 

Of 

THE 

INTEGRATION 

POINT 

IN 

THE 

GAMMA  DIRECTION. 

• CALL  CftLM 
•CALL  COLAS 
•CALL  COLA* 

•CALL  COLAS 
•CALL  CBLA7 
•CALL  CNOVA 
•CALL  CSLAN* 

C 

c 

DATA  T0L/5.0E-J/ 

C 

If (IP«R«T,CT.O)  SO  TO  300 

ifiRsr  • i 

EPOaSIBO/EL 

EPP  a 0,0 

CN2  a 0.0 

01002  a 1160**2 

CNl  a (0*5  - TNU)/EL 

CNS  a 1.0/EL 

TNUSQ  a TNU**2 

CN4  a (l.o  • T NU8S) *«2 

CN12  a (|.0  • TNU  ♦ TNUSQ) /CNR 

CNl 3 a (1,0  • 4,0*TNU  ♦ TNU83)/CN4 

CN«  a 0.75/ ((1.0  ♦ TNU) **2) 

CNS  a EL/ (1.0  • T NU**2) 

CN7  a EL*O.S/ (1*0  ♦ TNU) 

CNS  a 1,0/CN? 

LC  a 0 
LCMAX  a 100 
00  100  Lal,LMAX 
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AD-A035  644 
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unclassified 
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NL 


3cf3 


END 


3 -77 


noo 


l 


f 


IF  IMELF.fcfl.U  60  TO  70 
AlXX(l)  « 0.0 
ALTHL)  « 0,0 
ALXT(U  « 0.0 
iCl(UI)  • 0.0 
612 (L0  • 0.0 
8E3CL0  » 0.0 
T T MU  (L<>  ■ TNU 
70  EPBOUOaO.O 
100  KY(U  • 1 
C 

300  K$UM  a 0 

C 

IJ  a 

DO  3050  Kal.LBAR 
l a IJ  ♦ K 
HI  a ZF  IF) 

C DETERMINE  APPROPRIATE  REGION. 

KEY  a KY  CL) 

IF (KEY.GT.3)  60  TO  350 
60  TO  (400r*06'700)'  KEY 
J50  IF((K(V+l)/2.g8.KEY/2)  60  TO  600 
80  TO  700 

REGION  1.  ELASTIC  CURVE, 


400 


KSUM  a KSUM  ♦ 1 
IF (K8UH.GT . 1 ) 60  TO  450 

01  a CNb* (EXX  a TNU*£TT) 

02  a CNb* (ET  T ♦ TNU*EXX) 

03  a CN7*EXT 

04  a CNb»(XKXX  ♦ TNU*XKTT) 

05  a CNb* (XKTT  a TNU*XKXX) 
Ob  a CN7*XKXT 

81 A (H)  a 01 
82A(N)  a 02 
81A(M)  a 03 
84A (M)  8 04 
85A(N)  a 05 
64A(M)  a 06 
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ZZZ_^ __ _ — — ^ 


I 


490  SI  4 01  ♦ H|*04 

St  • oa  ♦ *1*09 

SS  • OS  ♦ Ht«St 

SISSO  I St • (61  • St)  ♦ 82**2  ♦ 9,0*69**2 

IF  (N«kP,t8,t>  60  TO  470 

CMOCkOtSISSO 

50  TO  IOOS 

17 0 IF(SI6SD,6B, 81801)  SO  TO  900 
BOSOM  • SISSO 
SS  TO  SOOO 

900  KTM  » KtV  ♦ I 

LBSBMtT  INTIRFOUTE  ON  518**  BAS  TO  CQARCCT  FOR  OVERSHOOT, 

SSSIS  4 SORT (SISSO) 

51  • SORT (COBQ(U) ) 


SI  • (S160  • Si) / (88616  • St) 

SI  • SIX  (L)  ♦ Si  * (01  • SXX(D) 

St  • STT(L)  ♦ St*(St  • STT(L)) 

89  4 SXT(t)  ♦ SI*(SS  • 9XT (t) ) 

SXSXXKL)  4 St 

SISTTI (t)  4 62 

BIBRTIU.)  4 69 
T|  4 CNS*(81  • TNU*6t) 

Tt  4 CN>*(6t  • TNU*6t) 

TI  4 CN9*6) 

EXX1M  ■ T1 
ETT1M  4 Tt 
IXTKU9  4 Tl 

IFS0  4 SSRT(CN|2* (Tl**t  ♦ T2**2)  • CN19*TI*T2  ♦ CN4*TI*»2) 

IPS(U)  • IPSO 

IFSOtkl  4 COSO 

IF  (JFIRST.tS.O)  JPIR8T  * 1 

SO  TO  9000 


RIS10NS  I AND  4,  DRASTIC  U0A01N6, 


HI  4 |XX  ♦ HI4RRXX  • SCI (l) 
N9  4 ITT  ♦ MUXKTT  • Slt(L) 
H4  4 (XT  ♦ HltXKKT  • SIS(L) 
CN*  4 TTNUU) 
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II  • 0 

blO  11  » II  ♦ 1 
CN22  ■ C>t2»t2 

EPBO  • SURTMU.Q  • CN2  ♦ CN(2)«(H2**2  ♦ Hl**2)  • 

1 (1.0  • «.0*CN2  ♦ CN(2)«H2*HS)/(1.0  • CMI2)»«2  ♦ 

2 0.75*M4**2/(1,0  ♦ CN2) **2) 

DELEP  • EPBD  • EPBO(L) 

EPP  • (EP«*DELCP  ♦ EU*CPBO(LO)/CPBO 
IP  IT  MU .6T.0.0)  C Hi  9 .5  • EPP*CN! 

IP  (AB8(CN2«TTNU(L)),LT. 0,000b)  60  TO  620 
IP  (II.6T.20)  60  TO  615 
TTNUU)  • CN2 
60  TO  610 

6)5  aRITE  (6.5500)  CN2,TTNU(L).TIME 

60  TO  4100 
620  CNR  « TTNU(L) 

IP  (EP60,LE«EPB0P)  BO  TO  650 
6i0  EP0CL)  • EPBO 

IP  (IPP.6T.EL.0R.EPP.LT.EP)  60  TO  4000 
81  ■ IPP/d.O  • CN2**2) 

62  • 0.5*IPP/(1.0  ♦ CNR) 

61  9 S1*(H2  ♦ QN2*MJ)  ♦ ALXX (L) 

62  • 81  * (NS  ♦ CN2*H2)  ♦ ALTT (L) 

6 J s S2*H4  ♦ ALXT(L) 

60  TO  1000 

c 

C 8ECONO  TEST  POR  UNLOAD IN6  IN  CITHER  RC6I0N  2 OR  4. 
650  01>X1A(M)*EXX1(L)  ♦ Hl*X4A(M) 

62>X2A(H)«CTTl (L)  ♦ Hl*X5A(H) 

QJM1A  (M) *EXT 1 (L)  ♦ HI *X6A (M) 

IP((P, (0*0.0)60  TO  660 
Pt«SXX(L)*8I6XRl(L)  ♦ ALXR(L) 

P2*8TT (L)*816TT1(L)  ♦ ALTT (L) 

PSMKT (L)*8I6XTl(L)  ♦ ALKT(L) 

60  TO  670 
660  P1«0.0 
P240.0 
PJaO.O 

670  E1401  • H2 
E2*02  • HI 
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11*91  • H4 

61*91*CN9*((1*TNU*(2) 

8iMi>CNI«((I*TNUt(U 

HMKNItll 

Al*9t*9t 

A2«9t*9t 

A 1*9 1*91 

»!B99*A|»(Al«A2)tAI«*2*),0*A)H| 

19  (9196D.6C,9160*,ANO.DELC9.8C,0.0)  60  TO  650 

C 

ftvuj*fttv#i 
TTNUU0*TNU 
E900(l) *91690 
911  (L9*Q1  ♦ 9(1(L) 

9E2U9*62  ♦ 9(2 (L) 

9(3(L0*91  ♦ 8EJCL) 

19(19,(8, 0,)80  TO  5000 
ALXX(69*9I 
AUTTCU0*9i 
ALXT(W0*91 

60  TO  5000 

c 

C 9(61 0M  5.  CLASTIC  UNL0AD1N6  • 8CL0A01N6. 

C 

700  Cl  • 9CI (L)  • EXX  • hl*XKXX 
E2  * 9(2(1)  • CTT  • HI *XK T T 
(5  • 9(1(L)  • (XT  • HlftXKXT 
Cl  • ALXX(L) 

C2  * ALTT(L) 

Cl  • AW«T(L) 

91  • Cl  - CN6*((1  ♦ TNU*(2) 

61  • Cl  • CN6*(M  ♦ TNU*(1) 

05  ■ C5  • CN7*C5 

A1  a 91  • Cl 

At  ■ 92  • Cl 

A1  • 95  - Cl 

91999  • A1*(A1  - A|)  * A2**2  ♦ 3.0*A1«*2 
19(91690, 6T, 91602)90  TO  900 
(980(l9*9169D 
60  TO  5000 


SIGMA 


] 90 


LlNUftlV  INTERFOLATE  ON  010MA  IAN  TO  CORRECT  FOR  OVEROMUOf* 
OUU  02  ■ SORT (EF90(L) ) 

60918  • SORT (91680) 

IF(02.6Tt0«O)  00  TO  646 
NC  • 0 

S20  01  • (90916  • 9160) / (00916  • 92) 

NC  ■ NC  ♦ 1 
IP(NC.6T.S)  60  TO  030 
OCll  • 91* (61  • 9XX (L) ) 

DEL*  • 01 • (62  • 9TT (L) ) 

OILS  • 91*(63  • 9XT (L) ) 

61  • 91  • 0EL1 

62  ■ 62  • DEL2 

63  ■ 63  • DELS 

A1  • 91  • ALXX(L) 

A2  « 92  • ALTT(L) 

A3  • 93  • ALXT(L) 

98016  • 9QRT(AI«(A1*A2)  * A2**2  ♦ 3,0*A3**2) 

IF (A 99 (96916*91 60) /91G0.6T , TOD  80  TO  029 
60  TO  635 

639  *RITE(fc,5700)  NC,K,I, J, KEY, S«0I6,§1 ,§2, TINE 
LC  • LC  ♦ 1 

IF  (LC.6T.LCMAX)  60  TO  4100 
935  CONTINUE 

DELI  • 61  • 9XX (L) 

DELE  * 62  • 8TTU) 

0EL3  > 63  • 8XT (L) 

T 1 » X1A(M)  ♦ CN3* (DELI  * TNU*0EL2)  * H1*X4A(H) 

T2  • X2A(M)  ♦ CN1* (0EL2  • TNU*DEL1)  ♦ Hl*X5A(N) 

TS  • X3A(M)  ♦ CN5*0EL3  ♦ H1*X*A(M) 

60  TO  990 

040  MRITEU'5200)  K,  I,  J,KET,  TIME, 02,90916 
T t « EXX  ♦ HltXKXX 

T2  * ETT  ♦ H1*XKTT 

T3  • EXT  ♦ HI *XKXT 

LC  ■ LC  ♦ 1 

IF  (LC.6T.LCHAX)  60  TO  4100 
000  EKXKLO  • T 1 
ETTKL)  • T 2 
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CBTlUft  * T9 
*»  • T 1 • 911(13 
H9  • T2  • IC2(L) 

N«  • TJ  • *11(13 

IMP  « 99RT<CNllt(Ml*«l  ♦ HI***)  • CN19*H|#H3  ♦ CMtNlMi) 
IPIOdP  • CMO 
8M(U)  • CMO 
918NX14L)  • 61 

9X6TT14L)  • 61 

*19*11(0  • 99 

KVCO  • Hit  ♦ 1 
90  TO  90*0 
C 

*000  9XX(L*  ■ 61 

*TT(tO  • *1 

9IT4L0  > 69 

>050  CONTI WOC 

K9UMA(M)  • KSUM 

c 

RETURN 

c 

C f**OR  RETURN, 

•000  HIT*  (*,5900)  EFF,«,I,J,TIMIfEF*D,CFB0P,EF90(l) 

•100  MtXTft  (6,5400) 

KERR  • 1 
1ITU1N 

5100  FORMAT  (I2H  IMMEDIATE  REL0A01N6  ,«l>, 3115,6) 

9900  FORMAT  (89M0EFF  19  OUT  OF  RAN9C,  CM  • 114,6/ 

1 119y  OCtl, A) 

9400  FORMAT  (UH090LUTION  X9  UN9TA9LC) 

5500  FORMAT  (19H  VALUE  OF  NU  WONT  CONVERSE, 2£l5. 6,  1SH  Tim,  9CC  * 

1 £15,6) 

5T00  FORMAT  (99H  CAN  NOT  TOTALLY  CORRECT  FOR  0VCR9H0QT/9IS*4E19.9) 
END 
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ftOECK  souve 

SUBROUTINE  SOLVE  ( A , N,  NDI M,  NOE  f , I f> , 01 T , B) 

C As  ORIGINAL'  MATRIX, 

c N a ACTUAL  DIMENSIONS  Of  A, 

C NO I M • OECLAREO  DIMENSION  OR  A IN  CALLING  PROGRAM, 

C NOIT  • DETC RMINCNT  COOS. 

C 0 * NOT  calculated, 

C 1 • CALCULATED. 

C IP  « INDEX  OP  K«TM  PIVOT  RON. 

C Oil  a DETERMINER!  OP  A, 

C S a RIGHT  HAND  SIDE  VECTOR, 

DIMENSION  A(NOIM,nO!mJ,IP(NOIM),B(NOIM) 

c ,<i.  iOiHi 

IPCNJal 

00  t Xai,N 

IP  4X.E0.N)  GO  TO  S 

XPl«X*t 

Max 

00  1 laKPi,N 

IP  (A»S(A(l.K)),Gr,A8SlA(M,K)))  Mai 

1 CONTINUE  ,*)) 

IP(R)aM 

IP  (M.NE.XJ  IP(N)a.IPlN) 

TaA(M,K)  if  i, 

A4N,X}aA(K,K) 

A (X, X) aT 

IP  IT. 10. 0.0)  GO  TO  3 
00  2 laxPtiN  TO 

2 A(IfK)a»A(I.X)/T 
00  « iaXPtiN 
?»A4M#I1 

A4M. 2) a A (X  # J ) 

A(X,J)af 

IP  (T.EG.O.O)  60  TO  4 
DO  3 laXPJ.N 

3 A(I,J)aA<I,J)*A(I,,.)*l 

4 CONTINUE 

5 IP  (a4K«K),EQ,0,0)  GO  TO  IS 

■j  it  ( . «• ) ,eu, o,o)  *o  ro  is 
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* CONTINUE 

IF  (NOET«EQ«t)  60  TO  li 
DET4IF4N) 

00  to  IM,N 

to  oiioiT*A(ifn 

11  IF  (NaEO,t)  00  TO  14 
nmi»n*i 

00  11  KMfNMl 
KFlSft*! 

M1F4R) 

TM(Mi) 

•IN)IB(R) 

0(K)Rf 

00  12  I*RF1 f N 

12  SU)M(I)tAU»K)«T 
00  IS  K6ai,NMi 
RMt«M*KO 
RlRNltl 

•<K)«»4K)/A(K(K) 

T«»64K) 

00  tl  lit $ RNl 
IS  0<I)M4»*A<I,K)*T 

14  B4|)9i4|)/A(lf  t) 

00  TO  17 

15  »RITI4t,14) 

16  FORMAT  (29H09XN6ULAR  MATRIX  IN  8/R  SOLVE) 
H f 0 

17  RETURN 
ENO 
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INITIAL  DISTRIBUTION 


Hq  USAF/SAMI  1 
AFIS/INTA  1 
ASD/F.NFFA  1 
Itq  PACAF/UOO  1 
AUL  (AUL/LSE-70-239)  1 
DDC  12 
Ogden  ALC/MMWM  2 
Hq  USAFE/DOQ  1 
AFATL/DLOSL  2 
AFATL/DL  1 
TAWC/TRADOCLO  1 
USNWC/Code  3263  1 
USNWC/Code  318  5 
USNWC/Code  326  1 
BRL  (AMXRD-BUL)  2 
BRL  (DRXBR-TE)  1 
ASD/XRO  1 
AFWL/SATS  1 
Waterways  Experiment  Station  2 
AFATL/DLYV  10 
ADTC/SD3M  1 
Hq  TAC/DRA  1 
Hq  USAFE/DOQ  1 
Hq  PACAF/DOO  1 
AFSC/XRPA  1 
Hq  USAF/RDQ  1 
Defense  Intelligence  Agency 

DB-4C3  1 
Hq  AFSC/SDZA  1 
Hq  TAC/XPSY  1 
Hq  TAC/DRFA  1 
Hq  TAC/DRAR  1 
AFFDL/FES  1 
AFAL/RW  1 
USAFTFWC/OA  1 
USAFTFWC/TE  1 
Hq  SAC/DOOB  1 
Hq  SAC/NR I (STINFO  Library)  1 
Comdr,  Naval  Wpns  Ctr/Code  403  1 
Comdr,  Naval  Wpns  Ctr/Code  317  1 
AFSC  Liaison  Office/Code  143  2 
Ogden  LAC/MMWMP  2 
Hq  TAC/DRA  1 
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